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ABSTRACT 


Control  of  a  modem  submarine  is  a  multi-dimensional  problem  coupling  initial 
stability,  hydrodynamic  and  control  system  response.  The  loss  of  stability  at  moderate 
to  high  speeds  is  examined  using  a  nonlinear  Hopf  bifurcation  analysis.  Complete  linear 
state  feedback  is  used  for  demonstration  purposes  for  depth  control  at  level  attitude  and 
for  a  fixed  nominal  speed.  Control  time  constant,  nominal  and  actual  speeds,  metacentric 
height,  and  stem  to  bow  plane  ratio  are  used  as  the  main  bifurcation  parameters.  A 
complete  local  bifurcation  mapping  provides  a  systematic  method  for  evaluating  the 
bounds  of  controllability  for  control  system  design  parameters  for  a  submarine  with  a 
given  set  of  hydrodynamic  coefficients.  The  submarine  and  its  potential  design 
modifications  are  then  verified  with  a  nonlinear  simulation  program. 
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closed  loop  dynamics  matrix  for  the  linearized 
system 
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(or,  control  surface  coordination  gain) 
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time  constant 

vehicle  forward  speed 
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vehicle  weight 

state  variable  vector 
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I .  INTRODUCTION 


The  fundamental  goal  of  submarine  control  is  to  reach  and 
maintain  ordered  depth.  Any  design  that  does  not  meet  this 
goal,  in  the  face  of  the  inherent  complexities,  is  not  overly 
useful  as  a  practical  vessel.  Current  evaluation  schemes 
involve  extensive  model  testing  such  as  those  done  on  the 
DARPA  SUBOFF  Model  (DTRC  Model  5470)  [ref.  8].  This  is  an 
expensive  and  time  consuming  evaluation  method.  The  goal  is 
to  develop  an  analytic  method  to  determine  the  stability  of  a 
potential  design.  Much  work  has  been  done  on  depth  control 
and  modelling  of  submarines  in  the  vertical  plane  [ref's.  1, 
2,  &  3]  but  this  work  has  assumed  the  existence  of  a  stable 
platform. 

The  stability  of  a  design  will  have  a  significant  impact 
on  its  responsiveness.  A  vehicle  with  a  large  margin  to 
instability  will  not  be  very  responsive.  The  problem  becomes 
one  of  determining  how  close  to  the  margins  we  can  get  without 
a  total  loss  of  stability.  Nonlinear  dynamics  and  chaos 
provides  us  with  the  tools  for  solving  this  problem  [ref's  9, 
&  10]  .  By  expanding  the  scope  of  our  research  to  include 
nonlinear  terms  we  are  able  to  define  the  limits  of  stability 
and  therefore  the  margins.  At  the  Naval  Postgraduate  School 
there  has  been  extensive  work  on  defining  the  nature  of  the 
instabilities  that  are  encountered  in  the  control  of 
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submersibles  in  path  keeping  and  the  dive  plane  [ref's  4,  5, 
&  6]  . 

At  low  forward  velocities  a  submersible  using  stern  planes 
for  attitude  and  depth  control  can  experience  a  loss  of 
stability  in  the  form  of  stern  planes  reversal.  This  is  a 
pitchfork  bifurcation  that  can  be  predicted  [ref.  7]  and  can 
be  accounted  for  in  the  control  design.  The  purpose  of  this 
thesis  is  to  develop  a  program  for  finding  the  limits  of 
stability  for  a  submarine  at  moderate  and  high  speeds.  Once 
these  limits  are  mapped  then  the  nature  of  the  loss  of 
stability  must  then  be  determined.  For  this  we  use  a  Hopf 
bifurcation  analysis  along  with  a  nonlinear  simulation 
program.  Finally,  after  the  limits  are  determined  we  are  able 
to  define  the  control  system  design  parameters  and  evaluate 
the  controllability  of  the  design. 
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II.  PROBLEM  FORMULATION 


A.  EQUATIONS  OF  MOTION 

By  restricting  the  submersible's  motion  to  the  vertical, 
or  dive  plane,  the  motion  can  be  modeled  by  the  coupled 
nonlinear  equations  for  pitch  and  heave.  With  a  body  fixed 
coordinate  frame  at  the  vehicles  geometric  center,  we  can 
express  Newton's  equations  of  motion  as 


m(w  -  Uq  -  zGq 


2  _ 


XgG)  =  Z4q  +  ZJW  +  ZqUq  +  ZuUw 

- cj*°SEbu » (r-xi|3dK 

Dj  TAIL  \w-xqi 

(W-B)c os0  +  U2(Zbtg  +  Z^S*)  , 


(2.1) 


and 


Iyq  +  mzGwq  -  mxG{w  -  Uq)  =  M^q  +  Mjw  +  MqUq  +  MwUw 

f  NOSE  l  W  —  XCf)  ^  , 

+CDf  b(x)- i; - ZQ—xdx 

DJ  TAIL  \w  -  XCt 


TAIL  \w  -  x$  (2.2) 

~(xGW  -  XgB)  c os0  -  ( zGW  -  ZgB)  sin0 

+  U2  (Mbbs  *  Mbbb)  . 


j.or  simplicity  we  will  assume  that  the  drag  coefficient,  CD, 
is  constant  over  the  length  of  the  submersible.  This  does  not 
have  a  significant  effect  on  the  results. 


3 


The  pitch  rate  and  depth  rate  for  a  submersible  is  given 
by 


b  =  -C7sin0  +  wcos0  ,  (2.4) 

respectively.  In  equations  (2.3)  and  (2.4)  8  is  the  pitch 
angle  of  the  vehicle  as  measured  from  the  horizontal  plane. 
Figure  (2.1)  shows  the  geometry  and  definitions  for  most  of 
the  symbols.  We  assume  that  forward  velocity,  U,  is 
maintained  constant  during  maneuvers  by  an  appropriate  use  of 
the  propulsion  system. 

B.  CONTROL  LAW 

Equations  (2.1)  through  (2.4)  can  be  written  as  a  set  of 
nonlinear  differential  equations  in  the  form. 


w  =  axlUw  +  a12Uq  +  a13zessin0  +  b21U2bg  +  b12U26b 
+  dv(w,q)  +  c2  ( w,  q)  , 


(2.6) 


g  =  a21Uw  +  a22Ug  +  a^z^inO  +  b21U2b£  +  b22U2hb 
*  dq(w,q)  +  c2  ( w,  q)  , 


(2.7) 


±  =  -C7sin0  +  wcos0  , 


(2.8) 


where, 

Dv  =  (m-Z*)  (Iy-Mj  -(mx^ZJ  (mXg+Af^)  , 
axlDv  =  (Iy-MJ  Zw+  (mxG+Z4)  Mv  , 
a12Dv  =  (Iy-Mj  (m+Zq)  +  (mxG+Zl))  (Mq-mxG)  , 
a12Dv  =  -  (mxG+Z4)  W  , 
blxDv  =  (Jy-A^)Z8<+(mxG+Z4)W4j)  , 
b12Dv  =  Uy-^jZ^Mmx^Z*)^  , 
a21Dv  =  (m-Zi,)  Mu+  (mxQ+M*)  Zw  , 
a22Dv  =  (m-Z*)  ( Mg-mxG )  +  (mxG+Af*)  (m+Zg)  , 
a23£»v  =  -  (m-Z*)  V  , 
b21Dv  =  (m-Z*)  (mxQ+Mj,)  Zb'  , 

bzz^v  ~  (m-ZJ)  ■Mjfc+  {tnxG+M^)  Zib  , 
dw(w,  q)  Dv  =  ( Iy-Mj  Iw+  (mxg+zj  Iq  , 
dq  ( w,  q)  dv  =  (m-zj  J9+  (mxG+Mw)  Iw  , 
cx  ( w,  g)  Dv  =  ( ly-Mf)  mzGq2  -  (mxG+Z^)  mzGwq  , 
c2  ( w,q)Dv  =  -  (m-Z^)  mzGwq+  (mx^Mj  mzGq2  . 


(2.9) 


In  equations  (2.5)  through  (2.8),  the  vehicle  is  assumed  to  be 
neutrally  buoyant  (W  *  B)  ,  level  (Xq  =  xB)  ,  and  statically 
stable  (zG  >  zB)  .  The  terms  Iw  and  Iq  represent  the  cross  flow 
drag  integrals  in  equations  (2.1)  and  (2.2),  and  zGB  =  zG  -  zB 
is  the  metacentric  height.  We  can  assume  zB  to  be  zero, 
therefore  zGB  ■=  zG. 

The  depth  control  system  uses  the  linearized  form  of 
equations  (2.5)  through  (2.8)  with  the  linearization  at  the 
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nominal  conditions  of  level  attitude,  ordered  depth  and 
forward  speed.  This  gives  the  state  space  form  of  the  depth 
control  equations  as 


6  =  q  , 


(2.10) 


w  =  allU0tr+a12U0q+a13zafl+b1U026  ,  (2.11) 


q  =  a21U0w+a22UQq+a23zGBe+b2u02t>  ,  (2.12) 


i  =  -U06+w  ,  (2.13) 

where  U0  is  the  nominal  speed  for  gain  selection,  and  the 
control  inputs  are  defined  as 

K  -  5  . 

bb-b*+ttb  (2.14) 

b\  *5ll  +  a*>12  • 

b2  ~  •^21+®^522  ' 

where  o i  is  the  control  surface  coordination  gain.  This 
produces  a  linear  full  state  feedback  control  law  of 

6  =  k1Q+k2w+kiq^kiz  ,  (2.15) 

where  the  gains  k2,  k2,  k3,  and  k4  are  computed  to  give  the 
closed  loop  system,  equations  (2.10)  through  (2.13),  the 
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desired  dynamics.  Since  the  desired  characteristic  equation 
has  the  general  form, 


X4+a3A,3+a2A,2+a1A.+a0  =  0  ,  (2.16) 

the  controller  gains  are  computed  by  equating  coefficients  of 
the  actual  and  desired  characteristic  equations, 

b\U<?  k2  +^2  ^02^3  =  _®3_  ^ll+<®22^  ^0  • 

b2U02k1+  (b2a12-bxa22)  UQ'k2+  (bxa2x-b2alx)  U^k3^bxU2kA 

=  ~a2~a23ZGB+  ^aiia22~a2iai2^  ' 

(2.17) 

{b2axx-bxa2x)  U02kx  +  (bxa23~b2ax3)  zGBU02k2 
+  {b2+bxa22-b2ax2)  C703/c4  =  ax+ (a13a21-a23axl)  zGBU0  , 

[  (b1a21~b2a11)  U04  +  (b1a23-b2a213)  zGBU0 2]  kt  -  cc0  . 

[ref.  7] 
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Figure  (2.1)  -  Geometrical  representations  of  the  basic 
definitions  for  the  equations  of  motion. 
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III.  STABILITY 


A.  BIFURCATION  ANALYSIS 

In  system  dynamics,  the  classical  definition  of  stability 
states,  that  the  real  part  of  all  the  eigenvalues  of  the 
system  must  be  negative.  Therefore,  our  initial 
investigations  into  the  stability  of  the  SUBOFF  Model  was  to 
find  those  eigenvalues  whose  real  parts  cross  the  imaginary 
axis.  We  used  the  bifurcation  analysis  program,  included  as 
Appendix  A,  to  calculate  the  eigenvalues  of  the  system. 

By  linearizing  the  equations  of  motion,  equations  (2.1 
through  2.4),  the  state  space  equations  of  the  dynamic  system 
can  be  written  in  the  form 


where 


u=-Kx 


(3.2) 


and  K  is  the  matrix  of  controller  gains,  as  calculated  by  pole 
placement  in  equation  (2.17).  The  eigenvalues  of  the  system 
are  found  by  solving 

(3  3) 

det\A-BK-sl\=0 .  v 
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In  the  bifurcation  analysis  program  a  pseudo- root  locus 
method  is  employed  where  the  time  constant,  Tc,  is  fixed.  The 
constant  Tc  fixes  to  placement  of  the  system  poles  at  a  given 
nominal  forward  speed  U0  and  then  the  model  speed,  U,  is 
varied  incrementally  with  the  system  eigenvalues  calculated  at 
each  speed  increment.  When  the  real  part  of  an  eigenvalue 
changes  sign  between  the  limits  of  a  speed  increment  a 
bisection  method  is  employed  to  find  the  speed  where  the  real 
part  of  the  eigenvalue  equals  0. 

For  each  point  where  the  real  part  of  an  eigenvalue 
crosses  the  imaginary  axis  the  associated  Tc  and  U  are  plotted 
on  a  bifurcation  map.  This  map  delineates  the  regions  of 
classical  stability  (all  eigenvalues  on  the  left  hand  half- 
plane)  from  the  regions  of  instability  (at  least  one 
eigenvalue  in  the  right  hand  half -plane) .  A  family  of 
bifurcation  maps  were  generated  by  varying  nominal  speed,  U0, 
initial  stability,  ZGB,  and  control  surface  gain,  a. 

B.  TYPICAL  RESULTS 

Figure  (3.1)  shows  a  typical  bifurcation  map  with  its  five 
distinct  regions.  Region  I  is  the  area  of  classical 
stability.  In  region  II  there  is  one  real  positive  eigenvalue 
which  is  indicative  of  a  pitchfork  bifurcation.  Pitchfork 
bifurcations  of  this  model  were  previously  examined  by  Reidel 
[ref.  7] .  Regions  III,  IV,  and  V  have  at  least  one  pair  of 
complex  conjugate  eigenvalues  with  a  positive  real  part.  This 
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would  indicate  that  there  should  be  an  unstable  oscillatory 
behavior  for  the  model . 

C.  SIMULATIONS 

An  extensive  set  of  simulations  were  run  to  verify  the 
bifurcation  map's  prediction  of  system  stability.  While  the 
results  of  simulations  showed  the  demarcation  between  the 
stable  and  unstable  regions,  the  simulations  found  that  the 
linear  bifurcation  analysis  failed  to  predict  the  method  of 
departure  from  stability.  Five  points  (a  through  e)  as  shown 
on  fig.  (3.1)  were  chosen  to  illustrate  the  model's  behavior 
in  the  regions  of  interest.  Table  (3.1)  lists  the  eigenvalues 
found  at  each  of  these  points  and  Table  (3.2)  shows  the 
eigenvalues  associated  with  the  exact  bifurcation  point  near 
points  (b  through  d) .  Note  that  the  eigenvalues  are  given  in 
dimensional  terms  while  all  other  information  in  the  tables 
are  non- dimensional ized. 

Point  a  is  in  the  region  of  stability  and  fig  (3.2)  shows 
a  rapid  convergence  to  nominal  stability.  Simulations  run  at 
points  b,  c,  and  d  are  shown  on  figures  (3.3),  (3.4),  and 
(3.5).  These  points  are  spaced  less  than  2.2  kts  apart  but 
they  have  a  huge  difference  in  their  bounded  oscillatory 
behavior.  Table  (3.1)  also  lists  the  measured  and  theoretical 
periods  for  these  three  points.  The  theoretical  period  is 
computed  by  dividing  2ir  by  the  imaginary  part  of  the 
eigenvalue  with  the  largest  real  part,  and  then  non- 
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dimensional ize  it  by  the  ratio  U/L,  where  U  =  speed  of  the 
vehicle  (in  feet  per  second)  at  each  point,  and  L  =  vehicle 
length  (13.9792  ft).  An  interesting  point  to  note  is  the 
behavior  at  point  d,  where  the  measured  period  is 
approximately  one  half  of  the  theoretical  period.  The 
dominant  period  at  80  nd  sec  is  associated  with  the  creation 
of  the  limit  cycle  at  the  bifurcation  point  while  the  155.4  nd 
sec  period  is  part  of  the  chaotic  response  found  at  point  d. 
Table  (3.2)  shows  the  eigenvalues,  measured,  and  theoretical 
periods  at  the  exact  bifurcation  points.  Finally  a  simulation 
run  at  point  is  shown  in  figure  (3.6)  and  demonstrates  an 
unbounded  departure  from  stability. 

It  is  evident  that  a  more  detailed  analysis  had  to  be 
performed,  and  that  a  Hopf  bifurcation  analysis  should  be 
used. 


TIME  CONSTANT  (non-dimensional) 


BIFURCATION  MAP  FOR  UO  =  9.0  FPS 


0  0.5  1  1.5  2  2.5  3 

VELOCITY  (non-dimensional) 


Figure  (3.1)  -  A  typical  bifurcation  map  showing  the  five 
distinct  regions . 
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Theta  (degrees) 


1 


COMPARATIVE  SIMULATION  PLOT 


Time  (nondimensional) 


Figure  (3.2)  -  An  example  of  the  stable  response  found  in 
region  I.  This  corresponds  to  Point  a  in  fig.  (3.1). 
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Theta  (degrees) 


COMPARATIVE  SIMULATION  PLOT 


Figure  (3.3)  -  This  figure  corresponds  to  Point  b  in  fig 
(3.1)  and  shows  the  oscillatory  behavior  in  region  III. 
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Theta  (degrees) 


COMPARATIVE  SIMULATION  PLOT 


Figure  (3.4)  -  This  simulation  corresponds  to  Point  c  in 
fig.  (3.1).  Note  the  distinct  difference  in  the  nature  of 
oscillations  between  this  point  and  Point  b. 
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COMPARATIVE  SIMULATION  PLOT 


Figure  (3.5)  -  This  simulation  corresponds  to  Point  d  in 
fig.  (3.1)  .  Note  the  magnitude  of  the  oscillations  associated 
with  region  IV. 
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COMPARATIVE  SIMULATION  PLOT 


Time  (nondimensionel) 


Figure  (3.6)  -  The  simulation  run  at  Point  e  in  region  V 
of  fig.  (3.1)  shows  the  unbounded  behavior  found  in  this 
region. 
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Table  (3.1)  -  A  listing  of  the  points  in  fig.  (3.1)  and 
their  associated  eigenvalues,  and  periods. 
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I  Points 

Eigenvalues 

U 

(fps) 

Measured  T 
(nd  sec) 

Theoretical  T 
(nd  sec) 

b 

-0.6070,  -0.0002 
0.0000  ±  0.4506 

0 . 6 1UQ 

6.11 

5.57 

c 

-0.5710,  -0.0054 
0.0000  ±  0.2884 

0.93U„ 

13.13 

CO 

rH 

• 

r> 

rH 

d 

-0.2675  ±  0.2063 
0.0000  ±  0.0558 

1 . 02Uq 

74.93 

74.67 

Table  (3.2)  -  Listing  of  the  eigenvalues  and  periods  of 
the  exact  bifurcation  points  associated  with  the  points  in 
Fig.  (3.1). 
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IV.  HOPP  BIFURCATION 


A.  INTRODUCTION 

By  definition  a  Hopf  bifurcation  occurs  when  a  pair  of 
complex  conjugate  eigenvalues  cross  into  the  right  hand  half- 
plane.  When  this  occurs  the  system  will  deviate  from  a  steady 
state  solution  in  an  oscillatory  manner.  This  deviation  can 
be  either  supercritical  or  subcritical. 

For  the  supercritical  case  shown  in  fig.  (4.1),  stable 
limit  cycles  form  after  straight  line  stability  is  lost. 
Assume  that  a  parameter,  w(t),  is  varying.  When  t  is  less 
than  tcrit  all  of  the  eigenvalues  of  the  system  are  located  in 
the  left  hand  half -plane  and  the  system  is  nominally  stable. 
At  t  equal  to  tcrit,  a  complex  conjugate  pair  moves  into  the 
right  hand  half -plane  and  forms  the  stable  limit  c/cle.  As 
the  distance  w(t)  increases  from  w(tcrit)  ,  where  the  difference 
D  =  w(t)  -  w(tcrit)  ,  the  amplitude  of  the  limit  cycle  will  also 
increase.  If  D  remains  small  then  the  system  will  remain  near 
the  nominal  steady  state  solution. 

Fig.  (4.2)  shows  the  subcritical  case  with  unstable  limit 
cycles  being  generated  prior  to  the  critical  point  being 
reached.  Thus,  as  w(t)  approaches  w(tcrit)  the  system  could 
deviate  from  the  nominal  steady  state  solution  and  converge  to 
a  large  amplitude  limit  cycle  before  the  nominal  system  loses 
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stability.  Hence  a  random  disturbance  can  cause  a  nominally 
stable  system,  one  with  all  of  the  eigenvalues  in  the  left 
hand  half -plane,  to  exhibit  oscillatory  behavior.  Once  w(t) 
equals  w(tcrit)  the  nominal  system  becomes  unstable  and  a 
discontinuous  increase  in  the  amplitude  of  oscillation  is 
seen.  This  is  because  there  are  no  nearby  stable  attractors 
for  the  system  to  converge  to. 

A  system  design  must  make  a  distinction  between  these  two 
types  of  bifurcations  because  of  the  disparate  nature  of  the 
losses  in  stability.  Thus  the  designer  cannot  rely  on  a 
linear  approximation  and  must  use  higher  order  approximations 
(3rd  order)  of  the  equations  of  motion  to  adequately  analyze 
his/her  system. 

B.  THIRD  ORDER  APPROXIMATIONS 

The  nonlinear  equations  of  motion  are  expanded  using  a 
third  order  Taylor  series  approximation  near  the  nominal 
steady  state,  x  =  [0] .  The  control  law  is  then  modelled  as, 

8'  «  6„etanh (-*£->  ,  (4-1> 
°sac 


where  5sat  is  the  saturation  angle  of  the  control  plane  input. 
Using  the  same  approximation  for  the  control  law  as  the 
equation  of  motion,  6'  becomes 
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1 


(4.2) 


b'  =  a- 


3«Mt2 


Therefore,  the  state  equations  can  now  be  written  as 

±  =  Ax  +  g{x)  ,  (4.3) 


where 


x=  [6  w  q  z]T  , 


(4.4) 


and  the  higher  order  terms  are, 

^1  =  0, 


(4.5a) 


gr2  =  bit72  6  3  (  6,w,g,z)  -  la13ZGB03 

6  (4.5b) 

+  cxlg2  +  c12wg  , 


g3  =  b2Uzb3  (0,  w,  q,  z)  -  \ a23ZGBQ 3 

6  (4.5c) 

+  c21wg  +  c22g2  , 


?4  =  +  7^3  •  (4.5d) 

A  6 
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The  term  $3  contains  the  third  order  terms  derives  from 
substituting  6  into  b' . 

Defining  T  as  the  matrix  of  eigenvectors  of  A  evaluated  at 
the  Hopf  bifurcation  point,  then  the  transformation, 

x  =  Tz  ,  (4.6) 


transforms  the  system  into  a  canonical  form 

&  =  T^ATz  +  T~1g(Tz)  .  (4.7) 


At  the  bifurcation  point 


T^AT 


0  -o)0  0  0 

Mq  0  0  0 

0  0  p  0 

0  0  0  q 


(4.8) 


with  o)0  >  0  and  p,q  <  0.  z3  and  z4  correspond  to  p  and  q  and 
are  asymptotically  stable.  Center  manifold  theory  states  that 
the  stable  coordinates  z3,  z4  can  be  expressed  as  polynomials 
in  the  critical  coordinates  z1#  z2  and  this  relationship  is  at 
least  second  order.  Therefore,  we  can  write, 

z^  =  a  1z12  +  a 2z^z2  +  a 2z22  ,  (4.9a) 


and 


*4  =  Pi*i2  +  P2ziz2  +  P3Z22  • 


(4.9b) 
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The  coefficients  ai  and  /3±  can  be  computed  as  follows:  We 
differentiate  equations  (4.9)  with  respect  to  z, 

£3  =  2a1z1£1  +  a2(z1z2  +  zxz2)  +  2a3z2£2  , 

£4  =  2PjZjij  +  P2  (£xZ2  *  ^P3^2^2  ' 

and  substitute  ±1  -  -w0z2  *  and  *2  *  “ozi  •  Therefore, 

£3  =  a2<i»0z12  +  (2a3  -  2 otj )  Ci)0z1z2  -  a2ti)0z22  ,  (4.10a) 


and 


£4  =  P2wo2i2  +  (2P3  *  2p  1)w0z1z2  -  P2«0z22  .  (4.10b) 


The  third  and  fourth  equations  of  (4.7)  are  written  as, 


3  =  |p  0|  z3 

4  lo  d  z4 


♦  [D]  , 


(4.11) 


where , 

[D]  =  T~1gz(Tz)  , 

and  g2(Tz)  contains  the  second  order  terms  of  g(Tz).  We 
substitute  equations  (4.9)  into  (4.11)  and  equate  coefficients 
with  (4.10)  .  In  this  way  we  get  a  linear  system  in  aif  and 
j8± .  From  this  we  can  write  the  two  dimensional  state  space 
equations  as 
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(4.12a) 


*1  =  "<*>0*2  +  -Til2!3  +  I12Z1^Z2  +  -^13^X^22  +  *14  *2* 
+  PUZX2  +  P12Zx22  +  Pi3Z22  , 


and 

*2  =  «0zl  +  r21zx3  +  r22z12z2  +  r23z1z22  ♦  r24z23 

(4.12b) 

+  P2l^l2  +  P22*l*2  +  P23^22  ' 

where  the  coefficients  rij  and  p^j  are  derived  from  equations 
(4.8)  . 

These  equations  are  only  valid  exactly  at  the  Hopf 
bifurcation  point.  For  speeds  U  in  a  region  near  the 
bifurcation  point  these  equations  become 

zx  =  a'ezj^  -  (<o0  +  w'e)  z2  +  F1(z1,z2)  ,  (4.13a) 


and 


t2  =  (»0  +  zx  +  a'ezj  +  F2(zx,z2)  ,  (4.13b) 


where:  a',  a>'  are  the  derivatives  with  respect  to  U  of  the 
real  and  imaginary  parts  of  the  critical  complex  conjugate 
pair  of  eigenvalues  evaluated  at  the  bifurcation  point:  e  is 
the  difference  in  U  from  the  critical  value;  and,  the 
nonlinear  functions  Fx  and  F2  are 
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(4.14a) 


z2  =  Fsin0  (4.15b) 

equations  (4.13)  become 

F  =  a'eR  +  Fx  (R, 0)  cos0  +  F2(F,0)  sin0  (4.16a) 

and 

F0  =  («  +  o'cJF  +  F2(F,0)  cos0  -  Fx (J?,0)  sin0  .  (4.16b) 

Equation  (4.16a)  then  yields 

F  =  a'cF  +  P(0)F3  +  0(0)  F2  .  (4.17) 

By  averaging  equation  (4.17)  over  one  cycle  we  can  obtain  an 
equation  with  constant  coefficients.  Defining 
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(4.18) 


k  =  -Lf2*P(e)dB  , 

2n  Jo 

and 

L  =  -Lf2*Q(9)  d&  ,  (4.19) 

2ltJo 

and  carrying  out  the  integration  we  obtain 

L  =  0  ,  (4.20) 

and 

K  =  (3ru  +  rl3  +  r22  +  3r24)  ,  (4.21) 

which  reduces  equation  (4.17)  to 

R  =  a'eR  +  iO?3  .  (4.22) 

The  existence  and  stability  of  the  limit  cycle  is 
determined  by  analyzing  the  equilibrium  points  of  the  averaged 
equation  (4.22),  which  correspond  to  periodic  solutions  in  zx 
z2  as  seen  in  the  coordinate  transformation  equations  (4.15) . 
From  equation  (4.22)  we  can  see  that  two  conditions  exist, 

1)  If  a’  >  0,  then: 

(a)  If  K  >  0  then  unstable  limit  cycles  coexist  with  the 
stable  equilibrium  for  e  <  0. 
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(b)  If  K  <  0  then  stable  limit  cycles  coexist  with  the 
unstable  equilibrium  for  e  >  0. 

2)  If  a'  <  0,  then: 

(a)  If  K  >  0  the  unstable  limit  cycles  coexist  with  the 
stable  equilibrium  for  6  >  0. 

(b)  If  K  <  0  then  stable  limit  cycles  coexist  with  the 
unstable  equilibrium  for  6  <0. 

From  this  criteria,  by  computing  the  nonlinear  coefficient 

K  we  can  use  it  to  distinguish  the  two  different  types  of  Hopf 

bifurcations : 

•  Supercritical  if  K  <  0; 

•  Subcritical  if  K  >  0. 

[ref.  6] 

C.  RESULTS 

A  typical  bifurcation  map  of  stability  for  the  SUBOFF 
Model  is  shown  in  fig.  (4.3) .  This  map  is  characterized  by 
the  Pitchfork  bifurcation  curve  (P)  and  the  three  Hopf 
bifurcation  curves  (HI,  H2,  and  H3)  .  The  nature  of  curve  P 
was  previously  analyzed  and  those  results  (Riedel,  1993)  are 
reconfirmed  in  this  study. 

Curve  HI  is  characterized  by  a  weak  supercritical  branch 
(a  -»  b)  at  low  nominal  speeds,  U0.  As  U0  increases  this 
branch  develops  a  weak  to  moderate  subcritical  behavior  with 
K  between  0  and  102.  The  second  branch  of  HI  (b  -»  c)  has  a 
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consistent  moderate  subcritical  behavior  with  K  on  the  order 
of  102. 

Cusp  (C)  marks  the  intersection  of  curve  H2  with  curve  H3 . 
The  cusp  is  highly  dependent  on  both  U0  and  initial  stability 

zgb- 

As  ZGB,  for  a  given  U0,  increases  curve  H2  (d  -*  e)  shifts 
from  a  very  weak  subcritical  nature  with  K  between  +10"2  and 
1  to  a  very  weak  supercritical  nature  with  K  between  -l  and  - 
10‘2.  With  a  lower  U0  and/or  higher  ZGB,  point  e  moves  down 
in  the  time  constant  and  may  not  intersect  curve  H3 . 

Curve  H3  (f  -*  g)  is  a  strong  supercritical  bifurcation 
with  K  values  between  -104  and  -106.  The  position  of  H3  is 
independent  of  U0,  initial  stability,  and  control  surface 
coordination  gain. 

Finally,  curve  H4  is  a  strong  subcritical  branch  with  K 
having  values  between  103  and  106.  Because  of  this  highly 
subcritical  behavior,  H4  can  dominate  and  obscure  the  stable 
region  at  speeds  greater  than  U/U0  =  1. 

Figure  (4.4)  plots  the  K  values  for  a  representative 
bifurcation  map  shown  in  fig  (3.1)  .  Note  the  predicted  super- 
and  subcritical  branches  associated  with  fig  (4.3)  .  Point  (a) 
is  inside  the  stable  region  I  and  the  numerical  simulations 
seen  in  fig.  (3.2)  converge  to  zero.  Point  (b)  is  located  in 
the  unstable  region,  immediately  after  a  supercritical 
bifurcation.  As  a  result,  small  amplitude  limit  cycle 
oscillations  have  developed.  The  same  is  true  as  we  move 
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towards  point  (c)  although  we  expect  a  family  of  unstable 
limit  cycles  around  this  point  as  a  result  of  the  subcritical 
bifurcation.  This  point  is  further  explored  in  the  next 
section.  As  we  approach  point  (d)  a  family  of  stable  limit 
cycles  is  generated  but  its  behavior  is  influenced  by  the 
previously  developed  unstable  limit  cycles.  The  real  part  of 
the  critical  pair  of  eigenvalues  is  becoming  positive  and 
relatively  large,  see  fig.  (3.1),  which  means  that  the 
amplitudes  of  these  stable  limit  cycles  are  expected  to  be 
significantly  higher,  a  result  which  is  confirmed  by  fig. 
(3.5).  The  imaginary  parts  of  the  critical  pair  of 
eigenvalues  are  also  changing  very  fast  in  this  region.  This 
means  that  the  period  of  these  limit  cycles  must  be  computed 
based  on  the  value  of  the  imaginary  part  right  at  the 
bifurcation  point,  rather  than  its  value  at  the  specific 
parameter  point.  This  was  observed  previously  in  tables  (3.1) 
and  (3.2).  Point  (e)  is  in  the  strongly  subcritical  region  V, 
thus  we  see  the  rapid  divergence  from  stability  as  seen  in 
fig.  (3.6). 

D.  SIMULATIONS 

The  response  of  the  submersible  was  simulated  using  an 
Adams -Bashf or th  integration  scheme  in  Fortran  program  coding 
and  is  included  in  Appendix  C.  The  contro7  law  (2.15)  and 
gains  (2.17)  discussed  in  Chapters  II  and  III  were  used.  Non- 
dimensional  ships  speed  and  control  system  time  constant. 
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nominal  speed,  initial  stability,  and  control  surface 
coordination  gain  were  the  input  values  to  tie  system.  A 
nominal  0.1  ft/sec  vertical  speed  was  used  as  the  external 
initial  disturbance. 

The  simulations  were  used  to  compare  the  Hopf  bifurcation 
data  in  two  ways: 

1)  By  confirming  sub-/  supercritical  behavior  predicted  by 
the  K  factor,  and; 

2)  By  comparing  the  predicted  to  simulated  period  of 
oscillations . 

In  fig.  (4.5)  we  have  plotted  the  stable  equilibria  and 
the  stable  and  unstable  limit  cycles  from  our  example  first 
used  in  fig.  (3.1).  Figure  (4.5)  clearly  shows  the  predicted 
sub-  and  supercritical  behavior  that  the  K  values  predicted. 
The  important  features  to  note  in  fig.  (4.5)  are: 

1.  Unstable  limit  cycles  found  with  the  subcritical  Hopf 
bifurcation; 

2.  The  divergence  of  tne  amplitudes  as  the  velocity  moves 
away  from  the  bifurcation  point  velocity,  and; 

3 .  The  rapid  divergence  of  the  right  most  bifurcation  and 
its  quick  and  abrupt  termination. 

To  see  why  the  abrupt  termination  occurs  we  must  look  at 
the  root  locus  plot  as  shown  in  fig.  (4.6)  .  The  parameter  in 
this  the  third  bifurcation  terminates  when  there  is  a  break  in 
point  and  the  complex  conjugate  pair,  which  had  given  the 
oscillatory  response,  becomes  two  positive  real  eigenvalues. 
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These  two  eigenvalues  now  make  the  system  completely  unstable. 
The  splitting  of  this  complex  conjugate  pair  into  two  real  and 
positive  eigenvalues  can  be  better  seen  in  fig  (4.7). 
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SUPERCRITICAL  HOPF  BIFURCATION 


Figure  (4.1)  -  An  example  of  a  supercritical  Hopf 

bifurcation  where  the  system  has  a  nominal  steady  state 
solution  until  point  C,  which  occurs  at  tcrit.  After  tcrit  the 
system  becomes  unstable  but  forms  a  stable  limit  cycle. 


34 


SUBCRITIC AL  HOPF  BIFURCATION’ 


Figure  (4.2)  -  This  figure  shows  a  subcritical  Hopf 

bifurcation  where  the  system  loses  stability  prior  to  reaching 
point  C,  which  again  occurs  at  tcric. 
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A  TYPICAL  BIFURCATION  MAP 


Figure  (4.3)  -  An  annotated  bifurcation  map  for  the  SUBOFF 
Model . 
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Figure  (4.4)  A  plot  of  the  K  values  associated  with  fig. 
(3.1) . 
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LIMIT  CYCLES  FOR  UO  =  9.0, Zg  =  0.4.  Alpha  =  0.0 


Velocity  (non-dimensional) 


Figure  (4.5)  -  This  shows  the  amplitude  response  for  a 
speed  range  encompassing  the  Hopf  bifurcation  points  that  are 
shown  in  fig.  (3.1) . 
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associated  with  fig.  (4.5). 


Real  Axis 


REAL  PART  OF  THE  EIGENVALUES  VS.  U 
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Figure  (4.7)  -  A  plot  of  the  real  parts  of  the  eigenvalues 
plotted  in  fig.  (4.6). 
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V.  APPLICATIONS 


A.  CONTROL  SYSTEM  PARAMETERS 

From  the  typical  bifurcation  maps  we  can  see  that  a  region 
of  stability  is  created  between  the  pitchfork  and  Hopf 
bifurcation  boundaries.  For  the  control  system  designer  the 
limits  or  parameters  must  be  defined  prior  to  starting  the 
design.  By  maximizing  the  region  of  stability  we  can  give  the 
designer  the  most  leeway  in  his/her  work.  There  are  three 
parameters  that  we  can  use  to  change  the  bifurcation  maps, 
nominal  speed,  initial  stability,  and  the  control  surface  gain 
coefficient. 

First  we  will  look  at  changing  the  nominal  speed.  In  fig. 
(5.1)  we  have  plotted  three  curves  for  nominal  speeds  of  3.0, 
9.0,  and  15.0  fps.  We  can  see  that  although  the  pitchfork 
line  moves  to  the  left  in  dimensional  speeds  this  line  remains 
nearly  constant  with  a  dimensional  stern  planes  reversal 
occurring  at  1.2  kts.  The  high  speed  Hopf  boundaries  (U/U0  > 
1.0)  move  apart  as  the  nominal  speed  increases.  The 
effectiveness  of  increasing  U0  is  limited  in  the  upper  branch 
by  the  fixed  position  of  the  H3  line  with  the  maximum 
practical  Tc  achieved  at  a  U0  =  9.0  fps .  In  the  lower  arm 
there  is  no  increase  in  the  stability  area  after  U0  =  9.0  fps 
therefore  any  increase  in  U0  is  pointless.  For  the  low  speed 
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Hopf  curves  (U/U0  <  1.0)  we  quickly  lose  our  margin  of 
stability  as  U0  increases  thus  necessitating  further  changes 
to  regain  the  lost  area  of  stability. 

The  next  parameter  change  we  examined  was  varying  the 
initial  stability.  Figure  (5.2)  shows  the  effect  of 
increasing  Zg  from  0.1  to  0.4  ft.  The  subcritical  H4  branch 
remains  constant  while  the  upper  high  speed  Hopf  branch  moves 
down  effectively  decreasing  the  area  of  stability.  The  low 
speed  Hopf  curves  move  up  to  increase  the  low  speed  area  of 
stability.  We  can  see  that  the  additional  loss  in  area  is  by 
the  movement  of  the  pitchfork  line  to  the  right.  At  a  Zg  = 
0.4  ft  stern  planes  reversal  occurs  at  a  dimensional  speed  of 
2.4  kts  which  is  well  within  the  currently  accepted  range  of 
1.0  to  3.0  kts  for  modern  submarines.  Therefore  in  this 
model  we  would  want  to  balance  the  initial  stability  to 
maximize  the  low  and  high  speed  areas. 

An  increase  in  the  control  surface  gain  coefficient,  a,  is 
shown  in  fig.  (5.3).  Note  that  the  low  and  high  speed  Hopf 
curves  all  move  up  in  time  constant.  While  the  low  speed  Hopf 
curves  give  a  large  increase  in  stability,  the  high  speed 
curves  move  up  proportionally  and  there  is  no  increase  in 
stability  area.  This  does  allow  the  designer  to  shift  the 
range  of  stable  time  constants  without  a  loss  of  high  speed 
stability.  The  pitchfork  line  will  move  to  the  left  until  it 
equals  zero  when  a  =  1.0. 
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In  order  to  examine  what  happens  at  the  extremes  of  the 
design  options  we  can  look  at  the  low  nominal  speed  (3  fps) 
bifurcation  maps.  Figure  (5.4)  shows  the  typical  bifurcation 
map  such  as  the  one  we  have  previously  discussed.  As  the 
metacentric  height  is  increased  there  are  significant  changes 
in  the  nature  of  the  bifurcation  curves.  In  fig.  (5.5)  we  see 
that  the  pitchfork  bifurcation  line  has  moved  significantly  to 
the  right  and  has  intersected  the  low  speed  Hopf  bifurcation 
curve.  This  intersection  along  with  the  merger  of  the  H2  and 
H4  curves  have  combined  to  reduce  the  region  of  stability  to 
a  negligible  portion  of  the  map. 

A  further  increase  in  the  metacentric  height  as  shown  in 
fig.  (5.6)  demonstrates  a  dramatic  change  in  the  nature  of  the 
stability  of  the  model.  The  low  speed  region  has  two 
hyperbolic- like  Hopf  bifurcation  curves  (the  upper  curve 
occurs  well  above  the  region  of  interest)  bounding  the  lower 
and  upper  limits  of  stability.  For  the  speeds  U/U0  >  1.0,  the 
pitchfork  bifurcation  line  now  intersects  the  H2  curve  and 
has  changed  from  a  supercritical  to  a  subcritical  pitchfork. 

Figure  (5.7)  is  an  example  of  the  effect  just  described 
and  how  it  can  occur  at  higher  nominal  speeds.  This  shows 
that  although  initial  stability  is  necessary  for  overall 
stability,  if  the  metacentric  height  becomes  too  large  it  can 
have  an  adverse  affect  on  the  performance  of  the  submarine. 


43 


B.  SUBMARINE  DESIGN  EVALUATION 

For  all  of  the  simulations  up  to  this  point  the  moment  and 
inertia  coefficients  of  the  control  surfaces  have  been  .5x  of 
those  listed  in  the  SUBOFF  Model  report.  By  using  the  reduce 
effect  control  surface  input  we  have  been  able  to  show  stable 
simulations  in  all  of  the  mapped  stable  regions  of  the 
bifurcation  maps. 

Linear  bifurcation  methods  fail  to  predict  the  change  in 
system  response  for  changes  in  control  surface  coefficients. 
The  bifurcation  maps  for  .lx  to  l.Ox  the  control  surface 


coefficients 

are  exactly 

the  same , 

in  other 

words 

the 

bifurcation 

points  are 

independent 

of  the 

size 

and 

effectiveness  of  the  control  surfaces.  Therefore  we  must 
examine  the  K  values  in  order  to  predict  the  response  of  the 
model . 

Figure  (5.8)  shows  the  change  in  stability  for  the  model 
with  and  increase  in  the  control  surface  coefficients.  The 
area  of  lost  stability  is  indicated  by  the  shaded  portion  of 
the  map.  This  loss  of  stability  is  caused  by  the  shift  of  the 
H2  curve  from  weak  to  moderately  supercritical  to  a  strongly 
subcritical  curve.  With  two  strongly  subcritical  curves  in 
the  high  speed  region  (U/U0  >  1.0)  the  possibility  of 
subcritical  capture  is  greatly  increased.  This  effect  is 
confirmed  by  running  an  extensive  set  of  simulations  and 
mapping  the  change  in  stable  to  unstable  response.  We  must 
note  that  this  instability  occurs  in  a  region  that  has  four 
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eigenvalues  with  negative  real  parts  where  linear  control 
system  design  would  not  predict  an  instability. 

With  this  method  a  design  can  be  tested  and  then  modified 
using  the  design  goals  as  determined  by  the  nonlinear  response 
in  the  simulations.  This  can  change  a  completely  unstable 
model  such  as  the  SUBOFF  Model  into  one  where  a  large  margin 
of  stability  and  control  system  latitude  is  design  into  it. 
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Time  Constant  (non-dimensional) 


COMPARATIVE  BIFURCATION  MAP:  INCREASING  Zg 


Figure  (5.2)  -  The  effects  of  changing  initial  stability, 
2g,  on  the  bifurcation  maps. 
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Time  Constant  (non-dimensional) 


COMPARATIVE  BIFURCATION  MAP:  INCREASING  Alpha 


Figure  (5.3)  -  The  effects  of  changing  the  control  surface 
gain  coefficient,  a,  on  the  bifurcation  maps. 


TIME  CONSTANT  (non-dimensional) 


STABILITY/BIFURCATION  MAP  FOR  UO  =  3  FPS 


Figure  (5.4)  -  A  low  nominal  speed  bifurcation  map. 
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TIME 


0  0.5  1  1.5  2 

VELOCITY  (non-dimensional) 


Figure  (5.5)  -  This  demonstrates  the  effects  of  an 

increase  in  the  metacentric  height  for  low  nominal  speeds. 
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BIFURCATION  MAP  FOR  UO  =  3.0  FPS 


N/Stable  Region 

O*  xaaa 


VELOCITY  (non-dimensional) 


Figure  (5.6)  -  As  the  metacentric  height  increases  it  has 
adverse  effect  on  the  stability  of  the  submarine. 


TIME  CONSTANT  (non-dimensional) 


BIFURCATION  MAP  FOR  UO  =  6.0  FPS 


Figure  (5.7)  -  An  example  of  increasing  the  metacentric 

height  for  a  higher  nominal  speed  model. 
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Figure  (5.8)  - 
moment  and  inertia 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

The  application  of  a  Hopf  bifurcation  analysis  to  a 
submarine  design  can  be  an  effective  tool  in  the  design 
evaluation  and  modification  phase.  These  methods  when  paired 
with  programs  that  generate  hydrodynamic  coefficients  for  a 
submarine  will  save  time,  effort,  and  money  by  reducing  the 
amount  of  model  testing  necessary  to  validate  a  design.  An 
effective  set  of  control  system  design  parameters  will  be 
generated  in  the  process  that  will  be  optimal  for  the  final 
design  of  a  submersible. 

Secondly,  this  system  will  give  the  limits  of  the  range  of 
metacentric  heights  that  will  maintain  stability  for  the  full 
range  of  speeds  of  the  design.  As  was  shown  changes  in  the 
metacentric  height  can  have  a  dramatic  affect  on  stability. 

Finally,  an  evaluation  of  the  need  for  bow  planes  or 
forward  control  surfaces  are  actually  to  maintain  control  of 
the  submarine.  If  the  forward  planes  can  be  eliminated  a 
potential  source  of  noise  would  also  be  eliminated  along  with 
simplifying  the  design  of  the  control  system. 
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B.  RECOMMENDATIONS 


First,  modify  the  programs  to  evaluate  the  performance  of 
the  submarine  for  the  effects  of  external  forces  such  as,  wave 
effects,  currents,  and  free  surface  effects. 

Second,  perform  a  systematic  study  of  the  bifurcations  at 
conditions  other  than  nominal  in  trim  conditions,  such  as  out 
of  neutral  buoyancy  (heavy/light) ,  or  out -of -trim  fore  and  aft 
conditions. 
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APPENDIX  A  -  BIFURCATION  ANALYSIS  PROGRAM 


PROGRAM  BIFUR1 . FOR 

BIFURCATION  ANALYSIS  FOR  THE  DARPA  SUBOFF  MODEL 

PARAMETERS  ARE:  TC  VS.  U 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  K1 ,  K2  ,  K3  ,  K4 ,  L , MQDOT, 

&  MWDOT , MQ , MW , MDS , MDB , MD , 

&  MASS, IY, WF,pi 

DIMENSION  A(4, 4) ,FV1(4) ,IV1(4) ,ZZZ(4,4) ,WR(4) ,WI(4) 

OPEN  ( 11 , FILE* ' BIF1 . RES ' , STATUS* ' NEW' ) 

OPEN  ( 12 , FILE* ' BIF2 . RES ' , STATUS* ' NEW' ) 

OPEN  ( 13, FILE* 'BIF3 .RES' , STATUS* ' NEW' ) 

WEIGHT*1556 .2363 


BUO 

*1556.2363 

L 

*13.9792 

IY 

-561.32 

G 

*32.2 

MASS 

-WEIGHT/G 

RHO 

=1.94 

XB 

*0.0 

ZB 

-0.0 

Pi 

*3.14159 

WRITE 

(*,*)  ' 

ENTER  MIN,  MAX,  AND  #  OF  INCREMENTS  IN  TIME 
CONSTANT' 

READ 

(*,*) 

TCMIN, TCMAX, ITC 

WRITE 

(*,*) 

'ENTER  MIN,  MAX,  AND  #  OF  INCREMENTS  IN 
VEHICLE  SPEED' 

READ 

(*,*) 

UMIN ,  UMAX ,  IU 

WRITE 

(*,*) 

'ENTER  NOMINAL  SPEED' 

READ 

(*,  *) 

UO 

WRITE 

(*,*) 

'ENTER  XG  AND  ZG' 

READ 

(*,*) 

XG,  ZG 

TCMIN 

*  0.1 

TCMAX 

*  9.0 

ITC 

*  200 

IU 

=  200 

UO 

H 

M 

to 

• 

o 

XG 

-  0.0 

WRITE 

(*,*) 

'ENTER  MIN,  MAX,  OF  VEHICLE  SPEED' 

READ 

(*,*) 

UMIN, UMAX 

WRITE 

(*,*) 

'ENTER  ZG' 

READ 

(*,*) 

ZG 
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c 

ZGB-ZG-ZB 
TCMIN-TCMIN*L/UO 
TCMAX-TCMAX* L/UO 
UMIN  =UMIN*UO 
UMAX  -UMAX*UO 
C 

ZQDOT--6.3300E-04*0.5*RHO*L**4 
ZWDOT--1.4529E-02*0.5*RHO*L**3 
ZQ  =  7. 54 5 OE -03*0. 5*RHO*L**3 
ZW  --1.3910E-02*0.5*RHO*L**2 
ZDS  -0.5M-5. 603  OE -03*0. 5*RH0*L**2) 

ZDB  -0.5* (-5.6030E-03*0.5*RHO*L**2) 

MQD0T--8 . 8000E- 04*0 . 5*RHO*L**5 
MWD0T--5 . 6100E-04*0 . 5*RHO*L**4 
MQ  «-3.7020E-03*0. 5*RHO*L**4 
MW  -  1 . 0324E- 02*0 . 5*RHO*L**3 
MDS  -0.5* (-2.4090E-03*0.5*RHO*L**3) 

MDB  -0 . 5*  {  2. 409 OE -03*0. 5*RHO*L**3 ) 

C 

ALPHA-- 0.0 
ZD- ZDS  + ALPHA* ZDB 
MD-MDS  +ALPHA*MDB 
C 

DV  = (MASS-ZWDOT) * (IY-MQDOT) 

&  - (MASS*XG+ZQDOT) * (MASS*XG+MWDOT) 

A11DV- (IY-MQDOT) *ZW+ (MASS*XG+ZQDOT) *MW 

A12DV-  (IY-MQDOT)  *  (MASS+ZQ)  +  (MASS*XG+ZQDOT)  *  (MQ-MASS*XG) 

A13DV-- (MASS*XG+ZQDOT) * WEIGHT 

BIDV  = (IY-MQDOT) *ZD+ <MASS*XG+ZQDOT) *MD 

A21DV-  (MASS-ZWDOT)  *MW+  (MASS*XG+MWDOT)  *ZW 

A22DV- (MASS-ZWDOT) * (MQ-MASS*XG) 

&  + (MASS*XG+MWDOT) * (MASS+ZQ) 

A23DV-- (MASS-ZWDOT) * WEIGHT 
B2DV  - (MASS-ZWDOT) *MD+(MASS*XG+MWDOT)*ZD 
C 

A11-A11DV/DV 
A12-A12DV/DV 
A13-A13DV/DV 
A21-A21DV/DV 
A22-A22DV/DV 
A23-A23DV/DV 
B1  -BIDV  /DV 
B2  -B2DV  /DV 
C 

EPS  -l.D-5 
ILMAX-1500 
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c 

DO  1  1*1 , ITC 

WRITE  (*,2001)  I, ITC 

TC-TCMIN+ (1-1) * ( TCMAX - TCMIN ) / (ITC-1) 

POLE-l.O/TC 
ALPHA3 *4.0* POLE 
ALPHA2-6 . 0*POLE**2 
ALPHA1-4 . 0*POLE**3 
ALPHAO*  POLE* *4 
K4-ALPHA0/ ( (B1*A21-B2*A11) *U0**4 
&  + (B1*A23-B2*A13) *ZGB*U0**2) 

A2M»B1*U0**2 

A3M-B2*U0**2 

A0M--  (A11+A22)  *U0-ALPHA3 
B1M*B2  *U0  *  *2 

B2M*  (B2*A12-B1*A22 ) *U0**3 
B3M*  (B1*A21-B2*A11)  *U0**3 

BOM* (A11*A22 -A21*A12 ) *U0**2 -A23*ZGB-ALPHA2 -B1*U0*U0*K4 
C1M= ( B2  * A1 1 - B1 * A2 1 ) *U0**3 
C2M* (B1*A23-B2*A13) *ZGB*U0**2 
COM*  (A13*A21-A23*A11) *ZGB*U0+ALPHA1 
&  - (B2+B1*A22-B2*A12) *K4*U0**3 

K2*C1M*B0M*A3M-B1M*C0M*A3M-C1M*B3M*A0M 
K2-K2/  (C1M*B2M*A3M-  B1M*C2M*A3M-  C1M*B3M*A2M) 

Kl- (C0M-C2M*K2) /C1M 
K3*  (A0M-A2M*K2)  /A3M 

DO  2  J*1,IU 

U-UMIN+  (J-l)  *  (UMAX-UMIN)  /  (IU-1) 

A(l, 1) =0 . 0D0 
A ( 1 , 2 ) *0 . 0D0 
A(l,  3) *1 . 0D0 
A ( 1 , 4 ) -0 . 0D0 
A(2, 1)  *A13*ZGB+B1*U*U*K1 
A(2 , 2)  *A11*U  +B1*U*U*K2 

A(2 , 3)  *A12*U  +B1*U*U*K3 

A (2 , 4)  *  B1*U*U*K4 

A  (3 , 1) «A23*ZGB+B2*U*U*K1 
A(3 , 2) *A21*U  +B2*U*U*K2 

A(3 , 3 ) =A22*U  +B2*U*U*K3 

A  (3 , 4 ) =  B2*U*U*K4 

A(4 , 1)  =-U 
A  (4 , 2 ) =1 . 0D0 
A(4, 3) =0 . 0D0 
A(4 , 4 )  *0 . 0D0 

1 

CALL  RG(4, 4, A,WR,WI, 0, ZZZ, IV1, FV1, IERR) 

CALL  DSTABL ( DEOS , WR , WI , FREQ ) 
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IF  (J.GT.l)  GO  TO  10 

DEOSOO-DEOS 

UOO  -U 

LL-0 

GO  TO  2 

DEOSNN-DEOS 

UNN  -U 

PR«DEOSNN*DEOSOO 

IF  (PR.GT.0.D0)  GO  TO  3 

LL-LL+1 

IF  (LL.GT.3)  STOP  1000 

IL-0 

UO-UOO 

UN-UNN 

DEOSO-DEOSOO 

DEOSN-DEOSNN 

UL-UO 

UR=UN 

DEOSL-DEOSO 
DEOSR-DEOSN 
U* (UL+UR) /2.D0 
A ( 1 , 1 ) =0 . 0D0 
A (1 , 2) -0 . 0D0 
A(1 , 3 ) =1 . 0D0 
A(l, 4) ®0 . 0D0 

A  (2 , 1 ) »A13*ZGB+B1*U*U*K1 
A (2 , 2 ) =A11*U  +B1*U*U*K2 

A(2,3)«A12*U  +B1*U*U*K3 

A  (2 , 4 )  =  B1*U*U*K4 

A(3 , 1) =A23*ZGB+B2*U*U*K1 
A  (3 ,2) »A21*U  +B2*U*U*K2 

A(3 , 3) =A22*U  +B2*U*U*K3 

A(3 , 4) =  B2*U*U*K4 

A(4, 1)  =-U 
A(4 , 2)  =1 . 0D0 
A  (4 , 3 ) =0 . 0D0 
A(4,4)*0.0D0 

CALL  RG (4 , 4 ,  A,  WR,  WI, 0 , ZZZ , IV1 , FV1 , IERR) 
CALL  DSTABL (DEOS , WR , WI , FREQ) 

DEOSM-DEOS 

UM-U 

PRL-DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF  (PRL.GT.O.DO)  GO  TO  5 

UO-UL 

UN*=UM 

DEOSO-DEOSL 

DEOSN=DEOSM 


IL-IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF-DABS (UL-UM) 

IF  (DIF.GT.EPS)  GO  TO  6 

U-UM 

GO  TO  4 

5  IF  (PRR.GT.0.D0)  STOP  3200 
UO-UM 
UN-UR 

DEOSO-DEOSM 

DEOSN-DEOSR 

IL-IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF-DABS (UM-UR) 

IF  (DIF.GT.EPS)  GO  TO  6 
U-UM 

4  LLL-10+LL 

WRITE  (LLL, 2002 )  U/U0,TC*U0/L 
3  UOO-UNN 

DEOSOO-DEOSNN 
2  CONTINUE 
1  CONTINUE 
C 

2001  FORMAT  (215) 

2002  FORMAT  (4 (3X, F10 . 6) ) 

END 

C 

SUBROUTINE  DSTABL (DEOS , WR , WI , OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR (4) ,  WI (4) 

DEOS= - 1 . OD+2  0 
DO  1  1-1,4 

IF  (WR(I) .LT.DEOS)  GO  TO  1 
DEOS-WR (I) 

IJ-I 

1  CONTINUE 
OMEGA- WI  (IJ) 

OMEGA-DABS (OMEGA) 

RETURN 

END 
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APPENDIX  B  -  HOPF  BIFURCATION  PROGRAM 


C 

C 

c 

c 


c 


c 


c 


PROGRAM  HOPF1A. FOR 

EVALUATION  OF  HOPF  BIFURCATION  FORMULAE  USING  THE  SUBOFF 
SUBMARINE  MODEL 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  L, IY, MASS , MQDOT, MWDOT, 

MQ , MW , MD , MDS , MDB ,  K1 ,  K2  ,  K3  ,  K4  , 

ALFA1 , ALFA2 , ALFA3 , 

BETA1 , BETA2 , BETA3 

DOUBLE  PRECISION  Mil ,M12 ,M13 ,M14 , M21 ,M22 ,M23 ,M24 , 

M31 , M32 , M33 , M34 ,M41 , M42 , M43 ,M44 , 

Nil , N12 , N13 , N14 ,  N21 , N22 , N23 , N24 , 

N31 , N32 , N33 , N34 , N41 , N42 , N43 , N44 , 

L21 , L22 , L23 , L24 , L31, L32 , L33 , L34 , 

L41 , L42 , L43 , L44 , L25 , L26 , L27, L35 , 

L36, L37, L21A, L22A, L23A, L24A, L31A, 
L32A, L33A, L34A 

DIMENSION  A (4,4) ,T(4,4) ,TINV(4,4) ,FV1(4) ,IV1(4) ,YYY(4,4) 
DIMENSION  WR ( 4 ) , WI ( 4 ) , TS AVE (4,4) , 

1  TLUD (4,4) , IVLUD ( 4 ) , SVLUD ( 4 ) 

DIMENSION  ASAVE (4,4) ,A1(4,4) ,A2 (4,4) 

OPEN  (10, FILE- 'HOPF. DAT' , STATUS- ' OLD ' ) 

OPEN  (20, FILE- 'HOPF. RES' , STATUS- ' NEW' ) 

WEIGHT-1556.2363 
IY  =  561.32 
L  -  13.9792 

RHO  -  1.94 

G  =  32.2 

XG  =  0.0 

ZB  =  0.0 

MASS  -WEIGHT/G 
BOY  -WEIGHT 

ZQDOT  =-6.3300E-04*0 . 5*RHO*L**4 
ZWDOT  =-1.4529E-02*0. 5*RHO*L**3 
ZQ  -  7. 545 0E -03*0. 5*RHO*L**3 
ZW  =-1.3910E-02*0. 5*RHO*L**2 
ZDS  =0.5* (-5.6030E-03*0.5*RHO*L**2) 

ZDB  =0.5* (-5.6030E-03*0. 5*RHO*L**2 ) 

MQDOT  =-8. 800 0E -04*0. 5*RHO*L**5 
MWDOT  =-5. 610 0E- 04*0. 5*RH0*L**4 
MQ  =-3. 702 0E -03*0. 5*RHO*L**4 
MW  =  1 . 0324E- 02*0 . 5*RHO*L**3 
MDS  =0.5*(-2. 409 0E -03*0. 5*RHO*L**3 ) 
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MDB 


0 . 5* (  2.4090E-03*0. 5*RHO*L**3 ) 


WRITE 

(*, 

1001) 

READ 

(*, 

*) 

ITOTAL 

WRITE 

(*, 

1002) 

READ 

(*, 

*) 

U0, ZG,DSAT 

WRITE 

(*, 

1004) 

READ 

<*, 

*) 

ZG 

WRITE 

(*, 

1003) 

READ 

(*, 

*) 

ALPHA 

ZGB-ZG-ZB 
UO  -9.0 
DSAT  -  0.4 

ZD- ZDS  + ALPHA* ZDB 
MD-MDS  +ALPHA*MDB 

DETERMINE  [A]  AND  [B]  COEFFICIENTS 


DV  - (MASS-ZWDOT) * (IY-MQDOT) 

- (MASS*XG+ZQDOT) * (MASS*XG+MWDOT) 

A11DV- (IY-MQDOT) *ZW+ (MASS*XG+ZQDOT) *MW 

A12DV-  (IY-MQDOT)  *  (MASS+ZQ)  +  (MASS*XG+ZQDOT)  *  (MQ-MASS*XG) 

A13DV-- (MASS*XG+ZQDOT) * WEIGHT 

BIDV  = (IY-MQDOT) *ZD+ (MASS*XG+ZQDOT) *MD 

A21DV- (MASS-ZWDOT) *MW+ (MASS*XG+MWDOT) *ZW 

A22DV- (MASS-ZWDOT) * (MQ-MASS*XG) 

+ (MASS*XG+MWDOT) * (MASS+ZQ) 

A23DV-- (MASS-ZWDOT) * WEIGHT 

B2DV  = (MASS-ZWDOT) *MD+ (MASS*XG+MWDOT) *ZD 


A11-A11DV/DV 
A12-A12DV/DV 
A13-A13DV/DV 
A21-A21DV/DV 
A22-A22DV/DV 
A23-A23DV/DV 
B1  -B1DV/DV 
B2  -B2DV/DV 

C11DV- (IY-MQDOT) *MASS*ZG 
C12DV-- (MASS*XG+ZQDOT) *MASS*ZG 
C21DV-- (MASS-ZWDOT) *MASS*ZG 
C22DV- (MASS*XG+MWDOT) *MASS*ZG 

C11-C11DV/DV 

C12-C12DV/DV 

C21-C21DV/DV 

C22-C22DV/DV 


o  n  o  non 


C 


DO  1  IT-1 , ITOTAL 

WRITE  (*,3001)  IT, ITOTAL 
READ  (10,*)  U, TC 
TC-TC*L/U0 
U-U*U0 

CALCULATE  THE  GAINS 

POLE-1 . 0/TC 
ALPHA3 -4.0* POLE 
ALPHA2 - 6 . 0 * POLE* *2 
ALPHA1-4 . 0* POLE* *3 
ALPHA0-  POLE* *4 
K4-ALPHA0/ ( (B1*A21 -B2*A11) *U0**4 
&  + (B1*A23-B2*A13) *ZGB*U0**2) 

A2M=B1*U0**2 

A3M=B2*U0**2 

A0M- - (A11+A22) *U0-ALPHA3 
B1M-B2*U0**2 

B2M=(B2*A12-B1*A22) *U0**3 
B3M= (B1*A21-B2*A11) *U0**3 

BOM-  (A11*A22  -A21*A12 )  *U0**2-A23*ZGB-ALPHA2-B1*U0*U0*K4 
C1M= (B2*A11-B1*A21) *U0**3 
C2M- (B1*A23 -B2*A13 ) *ZGB*U0**2 
COM-  (A13*A21-A23*A11)  *ZGB*U0+ALPHA1 
&  - (B2+B1*A22-B2*A12) *K4*U0**3 

K2=C1M*BOM*A3M-B1M*COM*A3M-C1M*B3M*AOM 
K2-K2/  (C1M*B2M*A3M-B1M*C2M*A3M-C1M*B3M*A2M) 

Kl= (C0M-C2M*K2) /C1M 
K3-  ( AOM  -  A2M*  K2  )  /A3M 


EVALUATE  NONLINEAR  RUDDER  EXPANSION  COEFFICIENTS 

DTTW-- (l.D0/(3.D0*DSAT**2) ) *3 .D0*K1*K1*K2 
DTWW=- (1 .DO/ (3 .D0*DSAT**2 ) ) *3 .D0*K1*K2*K2 
DTTQ=- (l.D0/(3.D0*DSAT**2) ) *3 .D0*K1*K1*K3 
DTQQ- - (1 .DO/ (3.D0*DSAT**2) ) *3 . D0*K1*K3*K3 
DTTZ-- (1 .DO/ (3 .D0*DSAT**2 ) ) *3 .D0*K1*K1*K4 
DTZZ-- (1 .DO/ (3 .D0*DSAT**2 ) ) *3 . D0*K1*K4*K4 
DWWQ-- (l.D0/(3.D0*DSAT**2) ) *3 .D0*K2*K2*K3 
DWQQ-- (1.D0/ (3.D0*DSAT**2) ) *3 .D0*K2*K3*K3 
DWWZ-- (1.D0/ (3.D0*DSAT**2) ) *3 .D0*K4*K2*K2 
DWZZ-- (1.D0/ (3.D0*DSAT**2) ) *3 .D0*K4*K4*K2 
DQQZ-- (1.D0/ (3.D0*DSAT**2) ) *3 .D0*K4*K3*K3 
DQZZ-- (1 .DO/ (3.D0*DSAT**2) ) *3 .D0*K4*K4*K3 
DTWQ-- (1.D0/ (3 .D0*DSAT**2 ) ) *6 .D0*K1*K2*K3 
DTWZ-- (l.D0/(3.D0*DSAT**2) } *6 .D0*K1*K4*K2 
DTQZ-- (l.D0/(3.D0*DSAT**2) ) *6 .D0*K1*K4*K3 
DWQZ-- (l.D0/(3.D0*DSAT**2) ) *6 .D0*K4*K2*K3 
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DTTT-- (1 .DO/ (3.D0*DSAT**2) ) *1 .D0*K1*K1*K1 
DWWW-- (l.DO/ (3 .D0*DSAT**2) ) *1 .D0*K2*K2*K2 
DQQQ-- (l.DO/ (3.D0*DSAT**2) ) *1 .D0*K3*K3*K3 
DZZZ-- (l.DO/ (3.D0*DSAT**2) ) *1 .DO*K4*K4*K4 

EVALUATE  TRANSFORMATION  MATRIX  OF  EIGENVECTORS 

A(1,1)-0.0D0 
A(1,2)-O.ODO 
A (1 , 3 ) -1 . 0D0 
A(1,4)«O.ODO 

A(2 , 1) -A13*ZGB+B1*U*U*K1 
A(2,2)-A11*U  +B1*U*U*K2 

A(2,3)«A12*U  +B1*U*U*K3 

A(2 , 4)  -  B1*U*U*K4 

A(3 , 1)  »>A2 3 * ZGB+B2 *U*U* K1 
A(3 , 2) »A21*U  +B2*U*U*K2 

A (3 , 3 ) =A22*U  +B2*U*U*K3 

A(3 , 4)  «=  B2*U*U*K4 

A  (4 , 1) * -U 
A(4 , 2)  *1 . 0D0 
A (4 , 3 ) =0 . 0D0 
A (4 , 4 ) =0 . ODO 
DO  11  1*1,4 
DO  12  J*1 , 4 

AS AVE ( I , J ) *A ( I , J ) 

12  CONTINUE 
11  CONTINUE 

CALL  RG (4,4 , A, WR, WI , 1 , YYY, IV1 , FV1 , IERR) 

CALL  DSOMEG ( IEV , WR , WI , OMEGA , CHECK ) 

OMEGAO = OMEGA 
DO  5  1*1,4 

T (I , 1) =  YYY ( I , IEV) 

T ( I , 2 ) *- YYY { I , IEV+1) 

5  CONTINUE 

IF  (IEV.EQ.l)  GO  TO  13 

IF  (IEV.EQ.2)  GO  TO  14 

IF  (IEV. EQ. 3 )  GO  TO  15 

STOP  3004 

14  DO  6  1*1,4 

T ( I , 3 ) * YYY (1,1) 

T (I, 4) =YYY (1,4) 

6  CONTINUE 
GO  TO  17 

15  DO  7  1-1,4 

T ( I , 3 ) = YYY (1,1) 

T (I , 4) =YYY (1,2) 

7  CONTINUE 
GO  TO  17 

13  DO  16  1-1,4 

T ( I , 3 ) -YYY (1,3) 
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T  ( 1 , 4 ) -YYY (1,4) 

16  CONTINUE 

17  CONTINUE 

CHECK  THE  EIGENVALUES 

WRITE (21 ,2003)  TC*U0/L,WI (1) ,WI (2) ,WI (3) ,WI (4) , 
&  WR (1) , WR (2) ,WR (3) ,WR (4) 

NORMALIZATION  OF  THE  CRITICAL  EIGENVECTOR 

INORM-1 

IF  (INORM.NE.O)  CALL  NORMAL (T) 

INVERT  TRANSFORMATION  MATRIX 

DO  2  1-1,4 
DO  3  J-1,4 

TINV (I, J) *0.0 
TSAVE(I, J)=T(I, J) 

3  CONTINUE 
2  CONTINUE 

CALL  DLUD (4,4, TS AVE , 4 , TLUD , IVLUD ) 

DO  4  1-1,4 

IF  ( IVLUD ( I ) . EQ . 0 )  STOP  3003 

4  CONTINUE 

CALL  DILU (4,4, TLUD , IVLUD , SVLUD ) 

DO  8  1*1,4 
DO  9  J=1 , 4 

TINV (I , J) =TLUD ( I , J) 

9  CONTINUE 

8  CONTINUE 

CHECK  Inv (T) *A*T 

IMULT*  1 

IF  (IMULT. EQ.l)  CALL  MULT (TINV, ASAVE,T,A2) 

IF  (IMULT. EQ.0)  STOP 
P=A2 (3,3) 

Q=A2 (4,4) 

WRITE (21,  *)  P,  Q 

DEFINITION  OF  Nij 

N11*TINV (1,1) 

N21-TINV (2,1) 

N31*TINV (3,1) 

N41*TINV (4,1) 

N12=TINV (1,2) 

N22-TINV  (2,2) 

N32-TINV (3,2) 
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N42«TINV(4,2) 

N13-TINV (1,3) 

N23-TINV(2,3) 

N33-TINV<3,3) 

N4  3  »TINV (4,3) 

N14-TINV (1,4) 

N24-TINV (2,4) 

N34-TINV (3,4) 

N44*TINV (4,4) 

C 

C  DEFINITION  OF  Mij 

C 

M11«T<1,1) 

M21-T (2 , 1) 

M31*T (3,1) 

M41«T (4,1) 

M12=T (1 , 2) 

M22=T (2 , 2) 

M32=T (3,2) 

M42*T (4,2) 

M13=T (1 , 3 ) 

M23=T (2 , 3 ) 

M3  3  =T (3,3) 

M43«=T (4,3) 

M14=T(1,4) 

M24=T (2,4) 

M34=T (3 , 4) 

M44=T (4,4) 

C 

C  DEFINITION  OF  Lij 

C 

L25=C11*M31*M31+C12*M21*M31 

C 

L26=2*C11*M31*M32+C12* (M21*M32+M22*M31) 

C 

L27=C11*M32*M32+C12*M22*M32 

C 

L35=C22*M31*M31+321*M21*M31 

C 

L36=2*C22*M31*M32+C21* (M21*M32+M22*M31) 

C 

L37=C22*M32*M32+C21*M22*M32 

C 

L21=  DTTW*M11*M11*M21  +  DTWW*M11*M21*M21  + 

&  DTTQ*M11*M11*M31  +  DTQQ*M11*M31*M31  + 

&  DTTZ*M11*M11*M41  +  DTZZ*M11*M41*M41  + 

&  DWWQ*M21*M21*M31  +  DWQQ*M21*M31*M31  + 

&  DWWZ  *M2 1 *M2 1 *M4 1  +  DWZZ*M21*M41*M41  + 

&  DQQZ  *M3 1 *M3 1 *M4 1  +  DQZZ*M31*M41*M41  + 

&  DTWQ*M11*M21*M31  +  DTWZ*M11*M21*M41  + 

&  DTQ  Z  *M1 1  *M4 1  *M3 1  +  DWQZ*M21*M41*M31  + 
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&  DTTT*M1 1  *M1 1  *M1 1  +  DWWW*M21*M21*M21  + 

&  DQQQ*M31*M31*M31  +  DZZZ*M41*M41*M41 

TTW-Mll*Mll*M22+2 . 0*M11*M12*M21 
TWW-M12*M21*M21+2 . 0*M11*M21*M22 
TTQ«Mll*Mll*M32+2 . 0*M11*M12*M31 
TQQ«M12*M31*M31+2 . 0*M11*M31*M32 
TTZ-M1 1  *M1 1  *M4 2+2 . 0*M11*M12*M41 
TZZ-M4 1  *M4 1  *M1 2+2 . 0*M11*M41*M42 
WWQ«M21*M21*M32+2 . 0*M31*M21*M22 
WQQ«M22*M31*M31+2 . 0*M31*M32*M21 
WWZ-M21*M21*M42+2 . 0*M41*M21*M22 
WZZ-M22*M41*M41+2.0*M41*M42*M21 
QQZ-M3 1  *M3 1  *M4  2 + 2 . 0  *M4 1  *M3 1  *M3  2 
QZZ-M32*M41*M41+2.0*M41*M42*M31 
TWQ-M1 1  *M2 1 *M3  2  +M3 1 *  (M11*M22+M12*M21) 
TWZ-M1 1  *M2 1  *M4 2  +M4 1  *  (M11*M22+M12*M21) 
TQZ-M1 1 *M4 1  *M3  2  +M3 1 *( Ml 1 *M4  2  +M1 2  *M4 1 ) 
WQZ-M2 1  *M4 1  *M3 2  +M3 1  *  (M2l*M42+M22*M41) 

TTT-3 . 0  *M1 1  *M1 1  *M1 2 
WWW-3 . 0*M21*M21*M22 
QQQ- 3 . 0  *M3 1  *M3 1  *M3  2 
ZZZ-3 . 0  *M4 1 *M4 1 *M4  2 

L22=DTTW*TTW+DTWW*TWW+DTTQ*TTQ+DTQQ*TQQ+ 

&  DTTZ  *  TTZ +DTZ  Z  *  TZ  Z  +DWWQ  *  WWQ +DWQQ  *  WQQ + 

&  DWWZ*WWZ+DWZZ*WZZ+DQQZ*QQZ+DQZZ*QZZ+ 

&  DTWQ*TWQ+DTWZ*TWZ+DTQZ*TQZ+DWQZ*WQZ+ 

&  DTTT*TTT+DWWW*WWW+DQQQ*QQQ+DZZZ*ZZZ 

TTW«M12*M12*M21+2 . 0*M11*M12*M22 
TWW-Mll*M22*M22+2 . 0*M12*M21*M22 
TTQ-M12*M12*M31+2 . 0*M11*M12*M32 
TQQ-M1 1  *M3 2 *M3 2+2 . 0*M12*M31*M32 
TTZ-M12 *M12 *M4 1+2 . 0*M11*M12*M42 
TZZ-M1 1 *M4 2  *M4 2+2 . 0*M12*M41*M42 
WWQ-M22*M22*M31+2 . 0*M21*M22*M32 
WQQ-M2 1  *M3 2  *M3 2 +2 . 0  *M2  2  *M3 1  *M3 2 
WWZ«M22*M22*M41+2 . 0*M21*M22*M42 
WZZ-M2 1 *M4  2  *M4  2+2 . 0*M22*M41*M42 
QQZ-M3 2*M32*M41+2 . 0*M31*M32*M42 
QZZ-M3 1 *M4  2  *M4  2+2 . 0*M32*M41*M42 
TWQ-M12  *M2  2  *M3 1 +M3  2  * (M11*M22+M12*M21) 
TWZ-M12  *M2 2  *M4 1 +M42  *  {Ml  1  *M2 2  +M12  +M2 1 ) 
TQZ-M1 2  *M4 2  *M3 1 +M3 2  *  (M11*M42+M12*M41) 
WQZ-M2 2 *M4 2  *M3 1 +M3 2 *  (M21*M42+M22*M41) 

TTT-3 .0*M11*M12*M12 
WWW-3 . 0*M21*M22*M22 
QQQ-3 . 0  *M3 1  *M3  2  *M3  2 
ZZZ-3 . 0*M41*M42*M42 
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L23«DTTW*TTW+DTWW*TWW+DTTQ*TTQ+DTQQ*TQQ+ 

&  DTTZ*TTZ+DTZZ*TZZ+DWWQ*WWQ+DWQQ*WQQ+ 

&  DWWZ*WWZ+DWZZ*WZZ+DQQZ*QQZ+DQZZ*QZZ+ 

&  DTWQ * TWQ +DTWZ *TWZ+DTQZ  * TQZ +DWQZ  * WQZ + 

&  DTTT*TTT+DWWW*WWW+DQQQ*QQQ+DZZZ*ZZZ 

C 

L24-  DTTW*M12*M12*M22  +  DTWW*M12*M22*M22  + 

&  DTTQ*M12*M12*M32  +  DTQQ*M12*M32*M32  + 

&  DTTZ*M12*M12*M42  +  DTZZ*M12*M42*M42  + 

&  DWWQ*M22*M22*M32  +  DWQQ*M22*M32*M32  + 

&  DWWZ*M22*M22*M42  +  DWZZ*M22*M42*M42  + 

&  DQQZ*M32*M32*M42  +  DQZZ*M32*M42*M42  + 

&  DTWQ*M12*M22*M32  +  DTWZ*M12*M22*M42  + 

&  DTQ  Z  *M1 2  *M4  2  *M3  2  +  DWQZ*M22*M42*M32  + 

&  DTTT*M12*M12*M12  +  DWWW*M22*M22*M22  + 

&  DQQQ*M32*M32*M32  +  DZZZ*M42*M42*M42 

C 

L31-L21 
L32=L22 
L33-L23 
L34-L24 

DEFINITION  OF  ALFA (I)  AND  BETA (I) 


D1  -N32*L25  +  N33*L35 
D2  =N32*L26  +  N33*L36 
D3  =N32*L27  +  N33*L37 
C 

Dll= - P 
D12=OMEGAO 
D21«=-2*OMEGAO 
D22=-P 
D23«2*OMEGAO 
D32--OMEGAO 
D33--P 
C 

ALFA2= (D2-D21*D1/D11-D23*D3/D33) / 

&  (D22 -D21*D12/D11-D23*D32/D33 ) 

ALFA1= (D1-D12*ALFA2 ) /Dll 
ALFA3- (D3-D32*ALFA2) /D33 
C 

D1  =N42*L25  +  N43*L35 
D2  «N42*L26  +  N43*L36 
D3  -N42*L27  +  N43*L37 
C 

D11--Q 

D12-OMEGAO 

D21«-2*OMEGAO 
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D22--Q 
D23-2*OMEGAO 
D32— OMEGAO 
D33--Q 
C 

BETA2- (D2 -D21*D1/D11-D23*D3/D33 ) / 

&  (D22 -D21*D12/D11-D23*D32/D33 ) 

BETA1- (D1-D12*BETA2) /Dll 
BETA3- (D3-D32*BETA2) /D33 
C 

L21A«2*C11*  (ALFA1*M31*M33+BETA1*M31*M34)  + 
C12*  (ALFA1*  (M2 1  *M3 3  +M2 3 *M3 1 )  + 

BETA1* (M21*M34+M24*M31) ) 

L2 2  A«2  *C1 1 * ( ALFA1 *M3  2  *M3  3  +BETA1 *M3  2  *M3  4 ) 
+2*C11* (ALFA3*M31*M33+BETA3*M31*M34) 
+C12* (ALFA1* (M22*M33+M32*M23) 

+BETA1* (M22*M34+M24*M32)  ) 

+012* (ALFA3* (M21*M33+M23*M31) 

+ BETAS*  (M21*M34+M24*M31)  ) 

L2 3A-2 *C1 1  *  (ALFA2 *M3 1  *M3 3  +BETA2 *M3 1  *M3 4 ) 
+2*C11* (ALFA3*M32*M33+BETA3*M32*M34) 
+C12*  (ALFA2*  (M21*M33+M23*M31) 

+BETA2* (M21*M34+M24*M31) ) 

+C12* (ALFA3* (M22*M33+M23*M32) 

+BETA3 *  (M22*M34+M24*M32)  ) 

L2  4  A®  2  *  Cl  1  *  (ALFA2  *M3  2  *M3  3  +BETA2  *M3  2  *M3  4 ) 
+C12* (ALFA2* (M22*M33+M23*M32) 
+BETA2*  (M22*M34+M24*M32)  ) 

L3 1A=2 *  C2  2  * ( ALFA1 *M3 1 *M3  3  +BETA1 *M3 1 *M3  4 ) 
+C21* (ALFA1* (M21*M33+M23*M31) 

+BETA1*  (M21*M34+M24*M31)  ) 


L32A=2 *C22 * (ALFA1 *M32  *M3  3  +BETA1 *M32  *M34 ) 
+2*C22* ( ALFA3  *M3 1 *M3  3  +BETA3  *M3 1 *M34 ) 
+C21* (ALFA1* (M22*M33+M32*M23) 

+BETA1* (M2  2  *M3  4  +M2  4  *M3  2 ) ) 

+C21* (ALFA3* (M21*M33+M23*M31) 

+BETA3  * (M21*M34+M24*M31) ) 

C 
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L33A-2*C22*  (ALFA2*M31*M33+BETA2*M31*M34) 
&  +2*C22* (ALFA3*M32*M33+BETA3*M32*M34 ) 

&  +C21* (ALFA2* (M21*M33+M23*M31) 

&  +BETA2* (M21*M34+M24*M31) ) 

&  +C21*  (ALFA3*  (M22*M33+M23*M32) 

&  +BETA3*  (M22*M34+M24*M32) ) 


L3  4A-2  *C2 2  *  (ALFA2  *M3  2  *M3  3  +BETA2  *M3  2  *M3  4 ) 

&  +C21* (ALFA2* (M22*M33+M23*M32) 

&  -*BETA2*  (M22*M34+M24*M32)  ) 

C 

L21=L21A+L21*B1*U*U-A13*ZGB*M11**3/6.D0 
L22«L22A+L22*B1*U*U-A13*ZGB*M12*M11**2/2.D0 
L23=L23A+L23*B1*U*U- Al3*ZGB*Mll*M12**2/2 .DO 
L24-L24A+L24*B1*U*U-A13*ZGB*M12**3/6.D0 
L3 1*L3 1A+L3 1  *B2  *U*U- A2  3  *  ZGB*M1 1 *  *  3 / 6 . DO 
L32«L32A+L32*B2*U*U-A23*ZGB*M12*M11**2/2.D0 
L33«L33A+L33*B2*U*U-A23*ZGB*M11*M12**2/2.D0 
L34«L34A+L34*B2*U*U-A23*ZGB*M12**3/6.D0 
L41«-0.5*M11*M11* (M21-U*Mll/3.0) 

L42- -Mil* (M12*M21+0 . 5*M11*M22 -  0 . 5*U*M11*M12 ) 
L43--M12* (Mll*M22+0 . 5*M12*M21- 0 . 5*U*M11*M12 ) 
L44=- 0 . 5*M12*M12* (M22-U*M12/3 . 0) 

Rll =N12 *L21+N13*L31+N14*L41 
R12=N12*L22+N13*L32+N14*L42 
R13=N12*L23+N13*L33+N14*L43 
R14=N12*L24+N13*L34+N14*L44 
R21=N22*L21+N23*L31+N24*L41 
R22=N22*L22+N23*L32+N24*L42 
R2 3  -N22  *L2  3  +N2  3  *L3  3  +N24  *L43 
R24=N22*L24+N23*L34+N24*L44 

EVALUATE  DALPHA  AND  DOMEGA 

UINC=0 . 001 
UR  =U+UINC 
UL  *=U-UINC 
U  =UR 

A(1,1)-0.0D0 
A (1 , 2 ) =0 . 0D0 
A (1 , 3 ) *1 . 0D0 
A(l,4) -0.0D0 
A  (2 , 1) =A13*ZGB+B1*U*U*K1 
A(2,2)»A11*U  +B1*U*U*K2 

A (2 , 3 ) =A12*U  +B1*U*U*K3 

A (2 , 4 )  *=  B1*U*U*K4 

A (3,1) =A23*ZGB+B2*U*U*K1 
A(3, 2) =A21*U  +B2*U*U*K2 
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A(3,3)-A22*U  +B2*U*U*K3 

A(3 , 4)  ■  B2*U*U*K4 

A(4,l)--U 
A(4,2)-1.0D0 
A(4 , 3) *0 . 0D0 
A(4 , 4) -0 . ODO 
C 

CALL  RG  (4 , 4 ,  A,  WR,  WI ,  0 ,  YYY,  IV1 ,  FV1,  IERR) 

CALL  DSTABL (DEOS,WR,WI, FREQ) 

ALPHR-DEOS 
OMEGRo FREQ 
C 

U-UL 

C 

A(l,  1)  »0 . ODO 
A(1,2)«0.0D0 
A(l, 3) *1 . ODO 
A(l,  4)  =*0 .  ODO 
A(2 , 1)  -A13*ZGB+B1*U*U*K1 
A  (2 , 2 ) =A11*U  +B1*U*U*K2 

A(2 , 3)  *=A12*U  +B1*U*U*K3 

A(2 , 4)  =  B1*U*U*K4 

A ( 3 , 1 ) =A2 3  *  ZGB+B2 *U*U*  K1 
A  (3 , 2 )  «=A21*U  +B2*U*U*K2 

A(3,3)«A22*U  +B2*U*U*K3 

A(3 , 4)  =  B2*U*U*K4 

A(4 , 1) =-U 
A (4, 2) “1 . ODO 
A  (4 , 3 ) =0 . ODO 
A(4,4)-0.0D0 
C 

CALL  RG ( 4 , 4 , A , WR , WI , 0 , YYY , IV1 , FV1 , IERR ) 

CALL  DSTABL (DEOS,WR,WI, FREQ) 

ALPHL-DEOS 

OMEGL*=FREQ 

C 

DALPHA= (ALPHR- ALPHL) / (UR-UL) 

DOMEGA- ( OMEGR - OMEGL ) / (UR-UL) 

EVALUATION  OF  HOPF  BIFURCATION  COEFFICIENTS 

COEF1-3 . 0*Rll+R13+R22+3 . 0*R24 
COEF2-3 . 0*R21+R23 -R12 -  3 . 0*R14 
AMU2  =-COEFl/ (8.0*DALPHA) 

BETA2=0 . 25*C0EF1 
C  TAU2  =- (COEF2-DOMEGA*COEFl/DALPHA) / (8 . 0*OMEGA0) 

PER  =2 . 0*3 . 1415927/OMEGAO 
PER  »PER*U/L 

WRITE  (20,2001)  TC*U0/L, COEF1 , DALPHA, 

&  PER , AMU2 , TAU2 , DOMEGA 

1  CONTINUE 
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STOP 

1001  FORMAT  ('  ENTER  NUMBER  OF  DATA  LINES') 

1002  FORMAT  ('  ENTER  U0,  ZG,  AND  DSAT' ) 

1003  FORMAT  {'  ENTER  BOW  PLANE  TO  STERN  PLANE  RATIO') 

1004  FORMAT  ('  ENTER  ZG' ) 

2001  FORMAT  (7E14.5) 

2002  FORMAT  (6E14.5) 

2003  FORMAT  (9F11.5) 

3001  FORMAT  (215) 

END 

SUBROUTINE  DSOMEG ( IJK , WR , WI , OMEGA , CHECK) 

IMPLICIT  DOUBLE  PRECISION  (A-HfO-Z) 

DIMENSION  WR (4 ) , WI (4 ) 

CHECK- - 1 . OD+2  5 
DO  1  1-1,4 

IF  (WR{I) .LT. CHECK)  GO  TO  1 
CHECK-WR { I ) 

IJ-I 

1  CONTINUE 

OMEGA-DABS (WI (IJ) ) 

IF  (WI(IJ) .GT.0.D0)  IJK-IJ 
IF  (WI(IJ) .LT.O.DO)  IJK-IJ- 1 
RETURN 
END 

SUBROUTINE  DSTABL (DEOS , WR , WI , OMEGA) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  WR { 4 ) , WI ( 4 ) 

DEOS- -1. OD+2 0 
DO  1  1-1,4 

IF  (WR(I) .LT.DEOS)  GO  TO  1 
DEOS-WR (I) 

IJ-I 

1  CONTINUE 
OMEGA- WI  (IJ) 

OMEGA-DABS (OMEGA) 

RETURN 

END 

SUBROUTINE  NORMAL (T) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  T ( 4 , 4 ) , TNOR (4,4) 

CFAC-T (1,1) **2+T (1 , 2 ) **2 

IF  (DABS (CFAC) .LE. (l.D-10) )  STOP  4001 

TNOR (1,1) =1 .DO 

TNOR (2,1) = (T(l,l) *T(2, 1) +T(2,2) *T(1,2) ) /CFAC 
TNOR (3,1)- (T (1, 1) *T (3,1) +T (3 , 2) *T (1, 2 ) ) /CFAC 
TNOR ( 4 , 1 ) « ( T ( 1 , 1 ) *T ( 4 , 1 ) +T ( 4 , 2 ) *T ( 1 , 2 ) ) /CFAC 
TNOR (1,2) *0 .DO 

TNOR (2 , 2 ) = (T (2 , 2 ) *T (1 , 1) -T (2 , 1) *T (1 , 2 ) ) /CFAC 
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TNOR(3,2)=(T(3,2)*T(l,l) -T (3 , 1) *T (1 , 2) ) /CFAC 
TN0R(4,2)=(T(4,2)*T(1,1) -T (4 , 1) *T (1 , 2 ) ) /CFAC 
DO  1  1-1,4 
DO  2  J-1,2 

T(I,J)=TNOR(I, J) 

2  CONTINUE 
1  CONTINUE 
RETURN 
END 


SUBROUTINE  MULT<TINV,A,T, A2) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  TINV (4,4) , A{4, 4) ,T(4, 4) , A1 (4, 4) , A2 (4, 4) 
DO  1  1-1,4 
DO  2  J-1,4 
A1 (I, J) =0 .DO 
A2  (I,  J)  *0  .DO 

2  CONTINUE 
1  CONTINUE 

DO  3  1-1,4 
DO  4  J-1,4 
DO  5  K-1,4 

A1 (I, J)-A(I,K) *T(K, J)+A1 (I, J) 

5  CONTINUE 
4  CONTINUE 

3  CONTINUE 
DO  6  1-1,4 

DO  7  J-1,4 
DO  8  K-1,4 

A2 (I, J) -TINV ( I , K) *A1 (K, J) +A2 (I , J) 

8  CONTINUE 

7  CONTINUE 

6  CONTINUE 
DO  11  1-1,4 

C  WRITE  (*,101)  (A(I, J) , J-1,4) 

11  CONTINUE 
DO  12  1=1,4 

C  WRITE  (*,101)  (T(I, J) , J-1,4) 

12  CONTINUE 
DO  10  1=1,4 

C  WRITE  (*,101)  (A2 (I, J) , J-1,4) 

10  CONTINUE 

C  WRITE  (*,101)  A2 (1, 1) 

RETURN 

101  FORMAT  (4E15.5) 

END 
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APPENDIX  C  -  ADAMS -BASHFORTH  SIMULATION  PROGRAM 

PROGRAM  ABSIM1 . FOR 

SUBOFF  SIMULATION  PROGRAM  USING  A  FOURTH  ORDER 
ADAMS - BASHFORTH  INTEGRATION  SCHEME 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

INTEGER  LDA,NA 
PARAMETER  (NA=4,LDA=4) 

DOUBLE  PRECISION  L,MW,MASS ,  IY,MQDOT,MWDOT,MQ,  ICP, LEN,  1 

*  MDS ,MDB, MDELTA, DELTA, UO , TII  3,TC,K1  K2 , 

*  K3,K4,F1{4)  ,F2(4)  (F3(4)  ,F%(4)  , 

*  A(LDA,NA)  ,  -SAT 
COMPLEX  EVAL(NA) 

DIMENSION  XL (25) ,BR(25) ,VEC1 (25) ,VEC2 (25) 

OPEN  (10 ,  FILE= '  SIM. DAT'  ,  STATUS  = '  OLD '  ) 

OPEN  (30, FILE- ' ANG . DAT ' , STATUS = ' OLD ' ) 

OPEN  (20, FILE=' SIM. RES' , STATUS* ' NEW' ) 

ENTER  SPEEDS  AND  TIME  DATA 

READ ( 10 , * )  UO , U, TC , TSIM, DELT, IPRNT, ALPHA 
U=U*U0 

GEOMETRIC  PROPERTIES  AND  HYDRODYNAMIC  COEFFICIENTS 

DSAT  =0 . 4000D0 

PI  =4 . 0*DATAN ( 1 . 0D0 ) 

RHO  =1 . 94D0 

CDZ  =0 . 5000DO 

L  =13 . 9792D0 

WEIGHT=1556 . 2363D0 
IY  =561 . 32D0 

MASS  =WEIGHT/32 . 2 
TC  =TC*L/U0 

ZQDOT=- 6 . 3300E- 04*0 . 5*RHO*L**4 
ZWDOT=-1.4529E-02*0.5*RHO*L**3 
ZQ  =  7.545 0E -03*0. 5*RHO*L**3 
ZW  =-l. 391 0E -02*0. 5*RHO*L**2 

ZDS  = (-5.6030E-03*0.5*RHO*L**2) 

ZDB  = (-5.6030E-03*0.5*RHO*L**2) 
MQDOT=-8.8000E-04*0.5*RHO*L**5 
MWDOT= -  5 . 6100E- 04*0 . 5*RHO*L**4 
MQ  =-3. 702 0E -03*0. 5*RHO*L**4 

MW  =  1.0324E-02*0 . 5*RHO*L**3 
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MDS  - (-2.4090E-03*0.5*RHO*L**3) 

MDB  -(  2.4090E-03*0. 5*RHO*L**3 ) 

C 

ZG  -0.4D0 

ZB  -0 . 0D0 

XG  -0 . 0D0 

XB  -0.0D0 

ZGB  -ZG-ZB 

ZD = ZDS + ALPHA* ZDB 
MD-MDS +ALPHA*MDB 

SET  INITIAL  CONDITIONS 

THETA=0 . 0D0 
READ (30,*)  THETA 
W  =0 . 1D0 

Q  =0 . 0D0 

Z  =0 . 0D0 

Z  =Z*L 

DETERMINE  [A]  AND  [B]  COEFFICIENTS 

DV  = (MASS-ZWDOT) * (IY-MQDOT) 

&  - (MASS*XG+ZQDOT) * (MASS*XG+MWDOT) 

AilDV= (IY-MQDOT) *ZW+ (MASS*XG+ZQDOT) *MW 

A12DV®  (IY-MQDOT)  *  (MASS+ZQ)  +  (MASS*XG+ZQDOT)  *  (MQ-MASS*XG) 
A13DV-- ( MAS S  *XG + ZQDOT ) * WEIGHT 
BIDV  = (IY-MQDOT) *ZD+ (MASS *XG+ ZQDOT) *MD 
A2 1DV= (MASS-ZWDOT) *MW+ (MASS*XG+MWDOT) *ZW 
A22DV* (MASS-ZWDOT) * (MQ-MASS*XG) 

&  + (MASS*XG+MWDOT) * (MASS+ZQ) 

A23DV= - (MASS-ZWDOT) * WEIGHT 
B2DV  = (MASS-ZWDOT) *MD+ (MASS*XG+MWDOT) *ZD 

A1 1 =A1 1DV/DV 
A12=A12DV/DV 
A13=A13DV/DV 
A21=A21DV/DV 
A22=A22DV/DV 
A23=A23DV/DV 
B1  =B1DV/DV 
B2  =B2DV/DV 

CALCULATE  THE  CONTROL  GAINS 

POLE-1 . 0/TC 
ALPHA3 -4.0* POLE 
ALPHA2 = 6 . 0 * POLE * * 2 
ALPHA1-4 . 0*POLE**3 
ALPHAO-  POLE* *4 
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K4-ALPHA0/ ( ( B1 *  A2 1 - B2  * A1 1 ) *U0**4 
&  +(B1*A23-B2*A13) *ZGB*U0**2) 

A2M=B1*U0**2 

A3M-B2*U0**2 

A0M= - (A11+A22) *U0-ALPHA3 
B1M«B2*U0**2 

B2M- (B2*A12 -B1*A22 ) *U0**3 
B3M- (B1*A21-B2*A11) *U0**3 

B0M= (A11*A22 -A21*A12 ) *U0**2 -A23*ZGB- ALPHA2 -Bl*U0*U0*K4 
C1M- { B2  *A1 1 - B1 *A2 1 ) *U0**3 
C2M= (B1*A23-B2*A13) *ZGB*U0**2 

COM- (A13*A21-A23*A11) *ZGB*U0+ALPHA1 
&  - (B2+B1*A22-B2*A12) *K4*U0**3 

K2  =C1M*B0M*A3M- B1M*C0M*A3M-  C1M*B3M* AOM 
K2-K2/  (C1M*B2M*A3M-B1M*C2M*A3M-C1M*B3M*A2M) 

Kl= ( COM- C2M*  K2 ) / C1M 
K3=  (AOM- A2M*K2  )  /A3M 

WRITE (*,*)  K1,K2,K3,K4 

DETERMINE  THE  EIGENVALUES 

A ( 1 , 1 ) =0 . 0D0 
A ( 1 , 2 ) =0 . 0D0 
A(l, 3) =1 . 0D0 
A ( 1 , 4 ) =0 . 0D0 

A(2 , 1) =A13*ZGB+B1*U*U*K1 
A (2 , 2 ) =A11*U  +B1*U*U*K2 

A  (2 , 3 ) =A12*U  +B1*U*U*K3 

A (2,4)=  B1*U*U*K4 

A(3 , 1) =A23  *  ZGB+B2  *U*U*K1 
A (3 , 2) =A21*U  +B2*U*U*K2 

A(3 , 3) =A22*U  +B2*U*U*K3 

A(3 , 4) =  B2*U*U*K4 

A(4, 1) =-U 
A(4 , 2) =1 . 0D0 
A(4 , 3 ) =0 . 0D0 
A (4 , 4) =0 . 0D0 

CALL  DEVLRG ( NA , A , LDA , EVAL ) 

CALL  DWRCRN  (' EVAL' , 1 , NA, EVAL, 1 , 0) 

DEFINE  THE  LENGTH  X  AND  BREADTH  B  TERMS  FOR  THE 
INTEGRATION 

XL (  1)=  0.0000D0 

XL (  2 ) =  0.1000D0 

XL (  3 ) =  0.2000D0 

XL (  4 ) =  0.3000D0 
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XL (  5)-  0.4000D0 

XL (  6)-  0.5000D0 

XL (  7)-  0.6000D0 

XL (  8)-  0.7000D0 

XL (  9)-  l.OOOODO 
XL (10)-  2.0000D0 

XL  (ID-  3 . 0000D0 
XL (12)-  4.0000D0 

XL (13)-  7.7143D0 

XL (14) *  10 . 0000D0 
XL (15)-  15 . 1429D0 
XL (16) »  16 . 0000D0 
XL (17) *  17 . 0000D0 
XL (18) «  18 . 0000D0 
XL (19)-  19 . 0000D0 
XL (20)-  20 . 0000D0 
XL (21)-  20 . 1000D0 
XL (22 ) -  20.2000D0 
XL (23)=  20 . 3000D0 
XL (24 ) =  20 . 4000D0 
XL (25)=  20 . 4167D0 

DO  102  N  =  1,25 

XL  (N)  =  (L/20 )  *XL  (N)  -  L/2 
102  CONTINUE 


BR (  1)=  0.000D0 
BR (  2 ) =  0.485D0 
BR (  3)-  0.658D0 
BR (  4 ) -  0.778D0 
BR (  5) =  0.871D0 
BR (  6 ) =  0.945D0 
BR (  7 ) =  1.010D0 
BR (  8) =  1.060D0 
BR (  9 ) =  1.180D0 
BR (10) =  1.410D0 
BR ( 11 ) =  1.570D0 
BR (12) =  1.660D0 
BR(13) -  1.670D0 
BR (14) =  1.670D0 
BR (15) =  1.670D0 
BR (16) -  1.630D0 
BR ( 17 ) -  1.370D0 
BR (18) -  0.919D0 
BR (19)-  0.448D0 
BR (20)-  0.195D0 
BR (21) -  0.188D0 
BR (22 ) -  0.168D0 
BR (23)=  0.132D0 
BR (24) -  0.053D0 
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BR (25)=  O.OOODO 
C 

PISIM  =TSIM/DELT 
ISIM  =PISIM 
C 

C  SIMULATION  BEGINS 

C 

DO  100  1=1, ISIM 
C 

C  CALCULATE  THE  DRAG  FORCE,  INTEGRATE  THE  DRAG  OVER  THE 

C  VEHICLE 

C 

DO  101  K=1 , 25 
UCF=W-XL(K) *Q 

VEC1 (K) =BR (K) *UCF*DABS (UCF) 

VEC2 (K) =BR (K) *UCF*DABS (UCF) *XL (K) 

101  CONTINUE 

CALL  TRAP (25,VEC1, XL, HEAVE) 

CALL  TRAP (25, VEC2 , XL , PITCH ) 

HEAVE =0 . 5*RHO*CDZ*HEAVE 
PITCH=0 . 5*RHO*CDZ*PITCH 
C 

C  COUPLING  EQUATIONS 

C 

DWDV= (IY-MQDOT) * HEAVE + (MASS*XG+ZQDOT) *PITCH 
DQDV= (MASS - ZWDOT) *PITCH+ (MASS*XG+MWDOT) * HEAVE 
C1DV= (IY-MQDOT) *MASS*ZG*Q**2 
Sc  -  (MASS*XG+ZQDOT)  *MASS*ZG*W*Q 

C2DV= - (MASS -ZWDOT) *MASS*ZG*W*Q 
&  + (MASS*XG+MWDOT) *MASS*ZG*Q**2 

C 

DW=DWDV/DV 

DQ=DQDV/DV 

C1=C1DV/DV 

C2-C2DV/DV 

C 

C  LIMITING  CONDITION 

C 

IF  (ABS (THETA)  .GT.  PI/4)  THEN 

WRITE  (*,*)  'BUBBLE  EXCEEDED  +/-  45,  BOAT  WAS  LOST' 
STOP 
ENDIF 
C 

C  CONTROL  EQUATIONS 

C 

DELTA  =  (K1*THETA+K2*W+K3*Q+K4*Z) 

DELTA  =  DSAT*DTANH (DELTA/DSAT) 

C  IF  (DELTA  .GT.  0.4)  DELTA=0 . 4 

C  IF  (DELTA  .LT.  -0.4)  DELTA=-0.4 

C 

THEDOT=Q 
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WDOT  «A11*U*W+A12*U*Q+A13*ZGB*DSIN (THETA) 
&  +B1*U*U*DELTA+DW+C1 

QDOT  -A21*U*W+A22*U*Q+A23*ZGB*DSIN (THETA) 
&  +B2*U*U*DELTA+DQ+C2 

ZDOT  --U*DSIN (THETA) +W*DCOS (THETA) 

IF  (I.GT.3)  GO  TO  150 

INITIAL  FIRST  ORDER  INTEGRATION 

THETA  -  THETA  +  DELT*THEDOT 

W  -  W  +  DELT*WDOT 

Q  -  Q  +  DELT*QDOT 

Z  -  Z  +  DELT*ZDOT 

FMD-THEDOT 
F2 (I) -WDOT 
F3 (I) -QDOT 
F4 (I) -ZDOT 

ADAMS -BASHFORTH  INTEGRATION 
150  FI (4) =Q 

F2 (4) -A11*U*W+A12*U*Q+A13*ZGB*DSIN (THETA) 

&  +B1*U*U*DELTA+DW+C1 

F3 (4) =A21*U*W+A22*U*Q+A23*ZGB*DSIN (THETA) 

&  +B2  *U*U*DELTA+DQ+C2 

F4 (4) — -U*DSIN (THETA) +W*DCOS (THETA) 

THETA- THETA+ (DELT/24.0) * (55.0*F1 (4) 

1  -59 . 0*F1 (3 ) +37 . 0*F1 (2) -9.0*F1 (1) ) 

W  -W+ (DELT/24.0) * (55. 0*F2 (4) 

1  -  59 . 0*F2 (3 ) +37 . 0*F2 (2) -9.0*F2 (1) ) 

Q  -Q+ (DELT/24.0) * (55. 0*F3 (4) 

1  -59 . 0*F3 (3) +37.0*F3 (2) -9.0*F3 (1) ) 

Z  -Z+ (DELT/24.0) * (55. 0*F4 (4) 

1  -59.0*F4 (3) +37.0*F4 (2) -9.0*F4(1) ) 

F1(1)-F1(2) 

FI (2 ) =F1 (3 ) 

FI (3) -FI (4) 

F2 (1) -F2 (2) 

F2 (2 ) -F2 (3) 

F2 (3 ) -F2 (4) 

F3 (1) -F3 (2) 

F3 (2 ) -F3 (3) 

F3  (3)  -F3  (4) 

F4 (1) -F4 (2) 

F4 (2) -F4 (3) 

F4 (3) -F4 (4) 

TIME-I*DELT 
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J-J+l 

IF  (J.NE.IPRNT)  GO  TO  100 
ALFA-W/U 

ALFA* (-DATAN (ALFA) )  * 180/PI 

WRITE  (20,991)  TIME*U0/L, Z/L, THETA*180/PI , 

&  DELTA, W,Q, ALFA 

WRITE  (*,991)  TIME*U0/L, Z/L, THETA* 180/ PI , 

&  DELTA, W,Q, ALFA 

J-0 
C 

100  CONTINUE 
STOP 

991  FORMAT  (2X, F8 . 2 , 6 (2X, F8 . 4) ) 

END 

C 

SUBROUTINE  TRAP (N, A, B, OUT) 

C 

C  NUMERICAL  INTEGRATION  ROUTINE  USING  THE  TRAPEZOIDAL  RULE 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DIMENSION  A ( 1 ) ,B(1) 

N1=N- 1 
OUT=  0 . 0 
DO  1  1*1, N1 

OUT1=0 . 5* (A (I) +A(I+1) ) * (B(I+1) -B(I) ) 

OUT  =OUT+OUTl 
1  CONTINUE 
RETURN 
END 
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APPENDIX  D  -  ROOT  LOCUS  PROGRAM 

%  THIS  MATLAB  PROGRAM  FINDS  THE  ROOT  LOCUS  FOR  A  GIVEN  UO,  Zg, 
%  Alpha,  AND  Tc. 

W-1556. 2363; Bl-1556. 2363 ; 

Iy-561.32; 
g«32 . 2 ; 
m=W/ g ; 
rho-1.94; 

L-13.9792; 
xg*0 ; zg« . 4 ; 
zgb* . 4 ; 

ndl=.5*rho*LA2; 
nd2= . 5*rho*LA3 ; 
nd3= . 5*rho*L^4 ; 
nd4=.5*rho*L^5; 

e  *  zeros (100,4) ; 

X  *  zeros (4,1) ; 

Y  »  zeros (4, 1) ; 

Z  *  zeros (4,1) ; 

Alpha=0 ; 

Zqdnd=-6 .33e-4;Zwdnd=-1.4529e-2;Zqnd*7.545e-3;Zwnd=-1.391e-2; 
Zds= .5*(-5.603e-3) ; Zdb=- . 5*5 . 603e-3 ; Zdltnd= (Zds+Alpha*Zdb) ; 
Mqdnda-8 . 8e-4 ;Mwdnd=-5 . 61e-4 ;Mqnd=-3 . 702e-3 ;Mwnd=l . 0324e-2 ; 
Mds=. 5* (-0.002409) ;Mdb= . 5*0 . 002409 ;Mdltnd= (Mds+Alpha*Mdb) ; 

Zqd=nd3*Zqdnd; Zwd=nd2*Zwdnd; Zq=nd2*Zqnd; Zw=ndl*Zwnd; 

Zdl t =ndl *  Zdl t nd ; 

Mqd=nd4  *Mqdnd ; Mwd=nd3  *Mwdnd ; Mq=nd3  *Mqnd ; Mw=nd2  *Mwnd ; 
Mdlt=nd2*Mdltnd; 

Dv= (m-Zwd) * (Iy-Mqd) - (m*xg+Zqd) * (m*xg+Mwd) ; 

allDva (iy-Mqd) *Zw+ (m*xg+Zqd) *Mw; 

al2Dv= (iy-Mqd) * (m+Zq) + (m*xg+Zqd) * (Mq-m*xg) ; 

al3Dv=- (m*xg+Zqd) *W; 

blDv= (Iy-Mqd) *Zdlt+ (m*xg+Zqd) *Mdlt; 

a21Dv= (m-Zwd) *Mw+ (m*xg+Mwd) *Zw; 

a22Dv»  (m-Zwd)  *  (Mq-m*xg)  +  (m*xg+Mwd)  *  (m+Zq)  ; 

a23Dv=- (m-Zwd) *W; 

b2Dv« (m-Zwd) *Mdlt+ (m*xg+Mwd) *Zdlt; 

all=allDv/Dv;al2-al2Dv/Dv;al3=al3Dv/Dv; 
a2 l=a2 lDv/Dv ; a2  2  »a22Dv/Dv ; a2  3  =a2  3Dv/Dv; 
bl=blDv/Dv;b2=b2Dv/Dv; 
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u  -  9; 

for  j  -  1:240 
U  -  j  /  2  0 ; 
si  (j)-  U; 

Fn-sqrt (UA2/ (zgb*g) ) ; 
tconst-4.7512; 
pole«-l/ (tconst*L/u) ; 

pl«[pole  pole-0.01  pole-0.02  pole-0.03]; 
a*[0  0  1  0;al3*zgb  all*u  al2*u  0;a23*zgb  a21*u  a22*u  0;... 
-u  1  0  0]  ; 

aU=  [0  0  1  0;al3*zgb  all*U  al2*U  0;a23*zgb  a21*U  a22*U  0;-U 
10  0]; 


b* [0;bl*uA2;b2*uA2;0] ; 
bU= [0  ;bl*U'*“2 ;b2*UA2 ; 0] ; 

Kl=place ( a , b , pi ) ; 

Y  =  eig (aU-bU*Kl) ; 
e  ( j  ,  :  )  =  Y'  ; 

end 

for  k  =  1:200 

U  -  8.9  +  k/1000 ; 
s2(k)  =  U; 

Fn=sqrt (U*2/ (zgb*g) ) ; 
tconst=4.7512; 
pole=-l/ (tconst*L/u) ; 

pl= [pole  pole-0.01  pole-0.02  pole-0.03]; 
a=[0  0  1  0;al3*zgb  all*u  al2*u  0;a23*zgb  a2l*u  a22*u  0;... 
-u  1  0  0]  ; 

aU= [0  0  1  0;al3*zgb  all*U  al2*U  0;... 
a23*zgb  a21*U  a22*U  0;-U  100]; 

b= [0;bl*uA2;b2*u^2;0] ; 
bU= [0;bl*UA2;b2*UA2;0] ; 

Kl=place (a,b,pl) ; 

Z  *  eig (aU-bU*Kl) ; 
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f (k,  :)  -  Z'  ; 
end 


axis ( [- .4  .1  - .8  .8] ) 
plot (real (e (1, 1) ) , imag (e (1, 1) ) , 'x' , . . . 
real (e(l,2) ) , imag (e (1, 2) ) , 'x' , . . . 
real (e(l,3) ) , imag (e( 1,3) ) , 'x' , . . . 
real (e (1,4) ) , imag (e (1,4) ) , 'x' ) 

hold 

plot (real (e ( : , 1) ) , imag (e ( : , 1) ) , ' . ' , . . . 
real (e ( : ,2) ) , imag (e ( : ,2) . 
real (e ( : , 3) ) , imag (e ( : , 3) . 
real (e(: ,4) ) ,imag(e(: ,  4)  )  , '  . ' ) 
plot ( real ( f ( : , 1 ) ) , imag ( f (:, 1 )),'.',.. . 
real ( f ( : , 2 ) ) , imag ( f (:, 2 )),'.',.. . 
real ( f ( : , 3 ) ) , imag ( f  (:, 3 )),'.',.. . 
real ( f ( : , 4 ) ) , imag ( f (:, 4 )),'.’ ) 
grid, title ('EIGENVALUES  FOR  SPEEDS  1  TO  100  FPS') 
%title(' DETAILED  LOW  SPEED  ROOT  LOCUS  PLOT') 
pause 
hold 
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